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ABSTRACT 


Direct-Sequence Spread Spectrum (DSSS) modulation is being advanced as the 
physical-layer basis for Seaweb undersea acoustic networking. DSSS meets the need for 
channel tolerance, transmission security, and multi-user access. This thesis investigates 
the performance of subspace-decomposition blind-equalization algorithms as alternatives 
to RAKE processing of DSSS signals. This approach is tailored for superior performance 
in time-dispersive and frequency-dispersive channels characteristic of ocean acoustic 
propagation. Transmitter and receiver structures are implemented in Matlab and 
evaluated with a statistics-based model of a doubly spread channel with additive noise. 
Receiver performance is examined using Monte Carlo simulation. Bit-error rates versus 
signal-to-noise ratio are presented for various multipath assumptions, noise assumptions, 
and receiver synchronization assumptions. 
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EXECUTIVE SUMMARY 


Underwater acoustic signalling is of increasing importance to the Navy since it 
extends network-centric warfare under the sea. In such applications small-scale channel 
effects like time variability and time spreading are significant when higher rate and 
coherent communications are desirable. Time spreading is caused by multipath 
propagation whereas time variability is the result of movement in the channel. Direct- 
Sequence Spread-Spectrum (DSSS) modulation is being advanced as the basis for 
undersea acoustic networking, because it meets the need for channel tolerance, 
transmission security, and multi-user access in the underwater environment. 

This thesis investigates the performance of subspace-decomposition blind- 
equalization algorithms as alternatives to RAKE processing of DSSS signals. These 
algorithms assume that the receiver has no a priori knowledge of the channel interference, 
but, since this interference spans a limited subspace in the eigenvalues of the 
autocorrelation matrix of the received signal, it can successfully be removed. This 
approach is tailored for superior performance in time-dispersive and time-spread 
channels. Transmitter and receiver structures are implemented in MATLAB and 
evaluated with a simulated statistics-based model of a doubly-spread channel with 
additive noise. Receiver performance is examined using Monte Carlo simulation. Bit- 
Error Rates (BER) versus Signal-to-Noise Ratio (SNR) are presented for various 
multipath assumptions, noise assumptions and receiver synchronization assumptions. 

The results suggest that the subspace decomposition method is suitable for 
underwater acoustic communications. A BER on the order of 10"^ is achievable within 
the desired SNR range. For a very severe channel, communications may be maintained at 
a lower bit-rate. 
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I. INTRODUCTION 


A. BACKGROUND 

The primary focus of this research is to improve a previously proposed Direct- 
Sequence Spread-Spectrum (DSSS) signaling scheme intended for use in Seaweb [1]. 
Seaweb is an experimental underwater acoustic network being developed by the Space 
and Naval Warfare Systems Center (SPAWARSYSCEN), San Diego which will extend 
network-centric warfare to the undersea environment [2], Seaweb is being developed to 
possess Low Probability of Intercept (LPI) and Low Probability of Detection (LPD) 
characteristics, high reliability, fault tolerance capability, and relative immunity to slow 
time-varying multipath. The previously proposed DSSS scheme combined a 
Differentially Encoded, spread In-Phase and Quadrature (IQ) components. Binary Phase 
Shift Keying (DS-IQ-BPSK) modulation which was proven very robust in the adverse 
underwater acoustic environment. The underwater acoustic channel effect has been 
simulated with an actual impulse response derived from the Signalex 2000 experiment 
(which means loss of generality) and by using the Bellhop propagation model. Bellhop is 
a static 2-D Gaussian ray model and is described in detail in [3]. Like most of the existing 
models for underwater acoustic channels. Bellhop is static in the sense that the impulse 
response only accounts for the multipath effect between transmitter and receiver, 
considering the surface to be an ideal flat reflecting boundary. However, the ocean 
surface is seldom stationary or smooth mainly because of the wind-induced surface 
waves. Part of the research in this thesis will focus on developing a more representative 
model that captures the effects of time variability. Based on this time-variable channel 
model, a better acquisition and equalization process will be implemented to further 
improve the DSSS scheme. 

B. SEAWEB REQUIREMENTS 

Under the current developments, Seaweb functions with node-to-node 
communication distances of three to five kilometers, using omni-directional transducers. 
The nodes are mounted on the bottom of the sea and no special geometry of the node 

distribution can be assumed because this varies according to the application. Moreover, 
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the network operates in shallow waters (50 to 200 m depth). The available operating 
frequency bandwidth is 9 to 14 kHz with provisions for a higher 15 to 22 kHz band. The 
acceptable bit-error rate (BER) is similar to the standards used in RF communication 
networks and is on the order of 10“^. Of interest to this thesis, Seaweb utility packets are 
fixed length sequences of 72 bits. The desired information transmission rate is 100 bits 
per second (bps). 

C. GOALS AND METHODOLOGY 

The first step in improving the signaling scheme is to use a more representative 
underwater acoustic model to account for the time variability and Doppler spread of the 
actual underwater channel. To introduce time variability, well known results and 
methodology, such as link-budget analysis and small-scale fading, are borrowed from the 
cellular communication community and adapted for the underwater environment. A static 
impulse response derived from a physics-based model (Bellhop) initiates the model. 
Subsequently, the static impulses are introduced with random x (or lag time) variability 
and, following a statistics-based approach, a new time-variable model is developed 
through independent Monte Carlo simulations and time averaging. Finally, the model 
allows the new time variable impulses to have a Gaussian width in order to simulate time 
spreading induced by scattering. This approach is not intended to produce a high-fidelity 
underwater acoustic model but rather an acceptable engineering model incorporating 
properties of the real underwater channel. 

The combination of the DS-IQ-BPSK and DSSS transmitting scheme was proven 
very robust as part of a previous thesis [1] as well as [4]. However, the RAKE receiver 
group based on the integrate-and-dump method was proven inadequate for coherent 
underwater communications. Therefore, a RAKE replacement is needed for the existing 
communication scheme in order to operate in the adverse underwater medium. This thesis 
focuses on developing an improved receiver structure based on advanced signal 
processing algorithms. The algorithms are implemented in MATLAB for various 
combinations of multipath, white and colored noise. Two blind-equalization algorithms 
are implemented based on whether synchronization between the transmitter and the 
receiver is achieved or not. The final step in the simulation is to test the developed 
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algorithms against the doubly-spread channel, using Monte Carlo simulation. All of the 
analysis takes place in the discrete time domain and the signals involved are considered 
to be impulses. Incorporating the algorithms developed here in the actual Seaweb modem 
will be the task of follow-on work. 

D. BENEFITS OF STUDY 

Underwater networks that employ wireless acoustic communication have a wide 
variety of applications in undersea warfare. Therefore, implementing reliable wireless 
underwater communications is of great importance to the Navy as it extends network¬ 
centric warfare under the sea. These wireless underwater acoustic applications include 
networked telemetry between sensors and base stations in littoral waters, submarine 
communications, control of minefields, and control of Unmanned Underwater Vehicles 
(UUVs). DSSS has been proven a durable and reliable scheme for underwater 
communication. However, the modem design is based on static models, which do not 
introduce time variability and Doppler spread like the real channel does. A time-variable 
model permits higher fidelity receiver design, making the existing scheme more robust. 

E. THESIS ORGANIZATION 

This thesis is organized into five remaining chapters. Chapter II discusses the 
large-scale sources of signal loss and attenuation in underwater communications and 
develops the associated link-budget analysis. Chapter III develops the small-scale fading 
effects of the underwater medium and analyzes the proposed statistics-based channel 
model. Chapter IV develops the theory for the DS-IQ-BPSK transmitter, and the 
simplified transmitter and receiver structures, emphasizing the signal processing 
algorithms implemented in the current communication scheme. Chapter V concerns the 
MATLAB implementation and performance measurements of the new communication 
scheme. Finally, Chapter VI reviews and summarizes the results and recommends follow- 
on work. 
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II. LARGE-SCALE FADING CHARACTERISTICS OF THE 


UNDERWATER ACOUSTIC CHANNEL 


This chapter discusses the large-scale fading effect, borrowing terminology from 
the digital communication community, of the underwater environment on the current 
communication scheme and uses link-budget analysis to predict the channel SNR and 
describe its performance. The underwater channel losses obey the one-way sonar 
equation and are formulated as an acoustic link-budget analysis. The overall losses and 
attenuation imposed on the transmitted signal by the medium are regarded to be the large- 
scale fading component. The severe fading of the pressure signal due to multipath (small- 
scale fading) in shallow water is addressed in the next chapter. 


A. ACOUSTIC LINK-BUDGETING 

The acoustic link-budget derives from the one way sonar equation: 

SNR = SL — TL — AN + transmitter receiver (2 -1) 


where all the quantities are in dB. The Signal-to-Noise Ratio (SNR) available at the 
receiver is defined as the ratio: 


SNR [dB] = 101ogio 


f-1 

sN, 


( 2 . 2 ) 


where S is the time average signal power and N is the time average noise power. Each of 
the remaining terms in Equation (2.1) is identified in the following paragraphs with a 
short description from [5]. 


I, Parameters Determined by the Medium 
a. Transmission Loss (TL) 

Transmission Loss (TL) is a metric used to describe the sum in dB 
(reference 1 pPa) of the signal loss from the transmitter to the receiver due to spreading 
and attenuation. Analytically this can be expressed as: 
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(2.3) 


TL = TL^+TL^ 

where TL^ is the component due to geometric spreading and TL^ is the component due to 
attenuation. 

Spreading is a geometrical effect that represents the weakening of the 
intensity of the sound field as the energy propagates into a larger volume. The signal is 
assumed to undergo a spherical spreading loss (inversely proportional tor^) to a range 
equal to the depth of the channel, gradually giving way to cylindrical spreading loss 
(inversely proportional to r). Spherical spreading loss is expressed as: 

Ji,,,,.,-, O'*) 

and the cylindrical spreading loss as 

TLg(cylindrical) = 10 ' logio ^ = 10- logjo (r) (2.5) 

where /j and are the time-average sound intensities in at the reference 

point (usually at 1 m from the source) and the receiver, respectively, and r is the range in 
[m]. 

Attenuation occurs as a result of absorption and scattering in the 
underwater channel. When sound propagates in the underwater medium part of the 
energy is transformed into heat and part is scattered in the ocean. Attenuation is used to 
describe the overall effect of leakage out of the sound channel by absorption and 
scattering as shown in Figure 1 on the next page. We observe that leakage dominates in 
the first region, ionic relaxation of two constituents of seawater dominates in the second 
and third, and the viscosity of seawater dominates the last region. 
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Region Region Region 

I II III IV 



Figure 1. Attenuation Coefficient a (/) [From Ref. 6.] 

The attenuation coefficient can be effectively calculated using Thorp’s 
expression [7] and each of the terms is respectively associated with the four regions 
described previously. 


«(/) = 3.3x10"'+0.11^^ + 44--^ + 2.75xlO"V" 

^ t rl A 1 r\r\ rl ^ 


( 2 . 6 ) 


! + /' 4100 + /' 

where / is the frequency in [kHz] anda(/) is the frequency-dependent absorption 
coefficient in [dB/km]. 

The TL^ component of TL can be calculated by: 


r4=a(/)T-io- 


(2.7) 


where r is the distance in [km] and a{f) is the frequency-dependent absorption 
coefficient in [dB/km]. 
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b. Ambient Noise (AN) 

Ambient Noise {AN) is the term describing the combined effect of oceanic 
turbulence noise, shipping noise, surface agitation noise and thermal noise at the receiver. 
AN is frequency dependent and different noise sources dominate the different frequency 
bands (in decades) [8]. A graphical representation can be seen in Figure 2 below. 



Figure 2. Ambient Noise (dB re IpPa) vs. Frequency (kHz) [After Ref. 8.] 

The first decade (1 Hz to 10 Hz) is dominated by oceanic turbulence 
whereas the second (10 Hz to 100 Hz) is dominated by shipping. For the next three 
decades (100 Hz to 100 kHz) wind-induced surface motion is the dominant feature and, 
for frequencies in excess of 100 kHz, thermal noise is dominant. In the frequency range 
used in Seaweb (9 to 14 kHz), we observe that the surface agitation noise is the dominant 
feature. 

Ambient noise has several kinds of variations: for example, diurnal, 
seasonal and geographical. Moreover, in shallow waters two intermittent noise factors 
must be considered: industrial noise (e.g., boats) and biological noise (e.g., snapping 
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shrimp). To incorporate these factors accurately in an analytical expression would require 
significant amount of in situ collected data for the deployment area of intent. 

The plot of the total ambient noise for different wind conditions in Figure 
3 below indicates that in the frequency range 9 to 14 kHz currently used in Seaweb, 
varying the wind speed from 0 to 20 m/s increases by 30 dB on the total ambient noise. 



Figure 3. Total Ambient Noise (dB re IpPa) vs. Frequency (kHz) for Various Wind 

Speeds [After Ref 8.] 


c. Large-Scale Fading Component ( 

The sum of all the above factors compose the overall signal degradation 
caused by the medium to the underwater acoustic channel. This sum is the large-scale 
fading component in wireless communication nomenclature and is measured 
indBre 1 pPa. Analytically: 

LSFC^,^„„^,=TL + AN. (2.8) 
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This intermediate result represents the total channel degradation that the signal must 
overcome for a given range. Figure 4 below is a plot of the channel SNR as a function of 
frequency for various ranges. 



The large-scale fading component increases with distance, which in turn means that the 
transmitted signal must compensate for more as the range increases. Furthermore, higher 
acoustic frequencies are impaired more significantly, implying that the available spectral 
bandwidth diminishes as a function of source-to-receiver range. 

2, Parameters Determined by the Equipment 
a. Directivity Index (DI) 

The directivity index (DI) is defined as the ratio of the intensity of a 
source at a specified direction divided by the intensity at the same reference point by an 
omni-directional source. In the case of Seaweb, the effect of DI is omitted since the 
transducers now in use are omni-directional. 
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b. Source Level (SL) 

The source level (SL) is defined [9] by the following expression: 


SL [dB re = 20■ log 


V2/2-Po 

V J 


(2.9) 


where i), is the peak acoustic pressure amplitude measured at a distance of 1 m from the 
source along its acoustic axis. In terms of power, the SL can be expressed as [5]: 

SL [t/5] = 171.5 + 10-log,oi^_,,, (2.10) 

where P is the power in [W] measured 1 m from the source. 


B. SUMMARY AND CONCLUSIONS 

The scope of the current chapter was to provide a general insight into the 
underwater acoustic environment for system design purposes. During the analysis, we 
observed that despite the similarities to the wireless networks, underwater 
communications are different from RF communications for several reasons. Firstly, 
underwater communications use acoustic pressure waves to propagate through the 
medium instead of electromagnetic waves. Secondly, the underwater medium strongly 
attenuates frequencies in excess of 30 kHz and therefore the bandwidth available for 
communications is very small compared to that of RF channels. Thirdly, the ocean noise 
is non-Gaussian. It is composed of discrete components, which can only be described 
analytically by “empirical” formulas like shipping noise, turbulence noise, surface noise 
and thermal noise, and other factors (e.g., industrial, biological, etc.) which are 
intermittent or can only be taken into account statistically if a large amount of data is 
available for a specific region of interest. 

Under the current link-budget analysis, a practical way of describing the large- 
scale fading effect in analogy to the cellular community was attempted. Several 
parameters have been unaccounted for resulting in a rough estimate of the link margin or 
SNR needed to communicate successfully underwater. If a more accurate estimate of the 
parameters in the sonar equation is necessary for a specific environment, one should 
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apply an underwater aeoustie propagation model to derive the losses. However, the link- 
budget will prove a handy tool for a general engineering approaeh to environmental 
losses and to ealeulate the required to establish communieation underwater. The 

interaetion with the medium boundaries and the assoeiated time variability and time 
spreading are eonsidered in the next ehapter. 
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III. SMALL-SCALE FADING CHARACTERISTICS OF THE 


UNDERWATER ACOUSTIC CHANNEL 


This chapter discusses the small-seale fading effects of the underwater 
environment. Two different charaeteristies of the underwater aeoustie medium are 
responsible for the small-seale fading: time spreading and time variability. Time 
spreading is caused by multipath propagation whereas time variability is the result of 
movement in the channel. As we see later in the ehapter time spreading introduees 
frequeney seleetivity of the channel while time variability introduees Doppler spread. It is 
not unco mm on for a real underwater channel to be simultaneously subjeet to time and 
Doppler spreading, a property that makes it a most unforgiving medium for 
eommunieation purposes. Figure 5 below depicts the time varying and spreading 
processes affeeting the signal. 



Figure 5. Processes Causing Small-Seale Fading Underwater 


An engineering model aeeounting for these effeets is required as an instrument for 
developing advaneed DSSS signaling algorithms. 
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A. 


SMALL-SCALE FADING CHARACTERISTICS 


I, Multipath 

The effect of multipath is that time-delayed echoes of the original signal arrive at 
the receiver at different times. The delay is caused by reflections from the sea surface and 
sea bottom. Multipath arrivals cause fluctuations in the received signal strength and the 
corresponding impulse response will have the general form: 

h{t) = '^w.{T)-5{t-T.) (3.1) 

1=1 

where rv, are the amplitudes associated with the i-th interaction and d is the Kronecker 
delta function used to denote the corresponding time-delayed version of the original 
signal. 

The effect of multipath in the transmitted signal is shown in Figure 6 below. The 
gross multipath delay is called and represents the duration of the overall impulse 
response. Values of in excess of ten milliseconds are not uncommon. For a transmitted 
signal with duration Ts, the received signal is no longer of duration T, but of duration Ts+ 
T 

± jfi , 


Ts 

Transmitted 


Channel 

Time 

Delay 

Spread 


<■ 


Ts + T„, 


> 


Received 


Figure 6. Effect of Time-Delay Spread on the Transmitted Signal 

Micropaths exist as a result of the various scatterers in the medium and on the 
reflective boundaries. The result of the reflection on a momentarily aligned boundary 
facet is that the unequal micropath lengths induce time spreading, randomize the phase 
and cause Rayleigh fading on the incident sound signal. Micropaths associated with sea- 
surface reflection are represented in Figure 7 on the next page. 
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Figure 7. Effect of Micropaths on Incident Signal 

The phenomenon as seen from the frequency-domain perspective is that the time 
delay spread introduces a bandwidth a l/T, (called coherence bandwidth) for which 

the channel passes all frequency components with equal gain. Within the coherence 
bandwidth, frequency components are well correlated and the channel transfer function is 
relatively flat. 

Based on the definition of the coherence bandwidth a channel may be 
characterized as either frequency selective or frequency non-selective based on the 
comparison of the signal bandwidth Ws to the coherence bandwidth of the channels^. 
Whenever the channel is called frequency selective and significant distortion 

occurs because the spectral components of the signal are attenuated in a frequency- 
dependent way. Since the above condition also implies thatT^ >r^, successive pulses 
overlap at the receiver causing inter-symbol interference (ISI). The spectral effect of a 
frequency selective channel on the signal is shown in Figure 8 on the next page. 
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Figure 8. Typical Frequency-Selective Fading When <W, 

The channel is characterized as frequency non-selective or flat when > W. In 
this case the various spectral components are equally affected by the channel. 
Furthermore, since this also implies that T^< f the channel will not induce ISI in the 
received signal. The spectral effect of a frequency non-selective channel is shown in 
Figure 9 below. 



2, Time Variability 

The time variability in the channel is introduced by the movement of the 
transmitter or the receiver or the channel itself. The sea-surface boundary is not static. It 
is the interaction with this boundary, combined with the motion of the source and 
receiver that accounts for most of the time variability of the real underwater channel. As a 
result, the impulse response of such a channel is not time invariant. In order to 
incorporate the time variability with time spreading. Equation (3.1) is rewritten as: 

= (3.2) 

(=1 
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where >v. (r) are the amplitudes associated with the i-th interaction and (t) is the 
associated time delay r , which is a function of time. 


The gross statistical measures descrihing the time variahility of the underwater 
channel is its coherence time and Doppler spread. By coherence time we mean the 
time duration during which the impulse response can be approximated as time invariant. 
In the frequency domain, the Doppler spread 5^ is the spectral broadening of the signal 

due to the channel variability. The two phenomena are related by « IjT ^. When both 
the transmitter and the receiver are fixed, the only time variability is caused by the 
surface wind induced waves and the corresponding surface motion. For telemetry systems 
the fluctuation of bandwidths is given by [10]: 


Bo=^-L- 


I I 4;r/ocos6>o 


(3.3) 


where w is the wind speed in [m/sec], /^=2/>v is the wave frequency in [Hz], 
= 0.005- is the wave height in [m], is the carrier (or center) frequency in [Hz], 
Oq is the incident grazing angle in degrees and c the speed of sound in [m/sec]. 


A Doppler spread channel is characterized according to the relation between the 
coherence time T^and the signal duration A channel is said to be fast fading when 

> Tp or equivalently B^>W where W is the bandwidth of the signal. Since 
W ^\/T^ = R^, where is the symbol rate, the channel impulse response changes 
rapidly during the time each symbol is transmitted. On the contrary, a channel is 
characterized as slow fading when or, equivalently, B^<W in which case the 

coherence time of the channel is greater than the transmitted symbol period and, 
therefore, less distortion is introduced as a result of the propagation in the medium. The 
effect of Doppler spread on a transmitted signal is shown in Figure 10 on the next page. 
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Figure 10. Effect of Doppler Spread on the Transmitted Signal 

3, Doubly Spread Channels 

A channel can display both time spreading and time variability simultaneously. 
Based on the relation of the product called the “spread factor”, a channel may be 

underspread if « 1 or overspread if • T,„ »1 .The most desired case is 

instinctively the under-spread condition for which a relatively high data rate can be 
achieved through the use of a coherent modulation method with coherent or differentially 
coherent detection. The definition of « 1 is always a point of discussion. In the 

case of Seaweb and the current communication scheme, a spread factor of less than 
10~^ is desirable. Given the multipath nature of the channel it is sometimes possible to 
resolve the multipaths effectively by reducing the spread factor. 


B. PROPOSED MODEL 

Existing underwater acoustic models are static and do not account for the small- 
scale factors. In such a model, the channel is idealized to impulses placed at lags r, with 

associated amplitudes iv^. and is always of the form of Equation (3.1). However, small- 
scale factors influence high-data-rate underwater communication schemes especially 
when a coherent scheme is required. The statistics-based model that follows is not 
intended to produce a high-fidelity underwater acoustics model but is aimed at providing 
a better simulation tool for analyzing communications performance. 

A static impulse response derived by Bellhop is shown in Figure 11 on the next 
page. The data were supplied courtesy of Dr. Paul Hursky of SAIC and are based on the 
exact bathythermograph conditions at which the Signalex 2000 sea trial was conducted 
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off the coast of San Diego. More details for the Signalex series experiments can he found 
in [2], Figure 11 helow displays the impulse response as a function of the source-receiver 
range. 
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Figure 11. Bellhop-Derived Series 0 to 6 km (Time vs. Range) 


The impulse response derived from Bellhop contained some very small 
amplitudes near the strong direct-path and multipath components. In order to have 
minimum overlapping when time variability is introduced, amplitudes smaller than 0.03 
were filtered out and only the dominant peaks were kept. The impulse response for range 
r = 5 km in magnitude vs. time axes is shown in Figure 12 on the next page: 
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Figure 12. Bellhop-Derived Impulse Response for r = 5 km 


We observe that the above impulse responses have exactly the form of Equation 
(3.1) and the lags r, are time invariable. In a time-varying channel, like the one presented 

in Figure 13 below, the lags r,. at which the multipath components occur are a function of 
time t. Therefore, if we snapshot the channel impulse response at different times t we 
observe different multipath components w, occurring at lags r, (0. 



Figure 13. Time Variability of a Cellular Communication Channel [After Ref. 11.] 
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In order to introduce time variability to an otherwise static impulse response, as 
the one received from Bellhop, we allowed each of the multipath components to vary 
randomly within an arbitrary range + 3 lags in x. The amplitudes associated with each 
lag were not altered at this stage and are the original amplitudes derived by Bellhop. The 
random phase shift caused by irresolvable micropaths was approximated by adding a 
uniform random phase factor from [0,27r] to each of the amplitudes. The process was 
simulated a thousand times in order to get a good statistical sample and the position of 
each lag was ensemble averaged in time t. This process is similar to collapsing the time t 
axes in Figure 13. The resulting impulse response is time varying, as each time lag x is 
now a function of time t and possesses a random phase shift. The new impulse response 
has the form: 

Kt, 0 = S (3-4) 

i=l 

where >v, (r) are the amplitudes associated with the i-th interaction, is the associated 
time delay x , which is a function of time, and the uniform phase shift. Furthermore, to 

better regulate the amount of time spread imposed by the real channel, the delta function 
was replaced by a Gaussian, except for the direct path component. To avoid rescaling the 
amplitudes, the factor l/crV^was omitted. The final form of the impulse response is: 

n 

- T ( 0 ) (3.5) 

/-I 

where g^{t) = Qyxp{-x^jlal) is a Gaussian shaped function with zero mean and cr3the 
desired time spread for multipath i. Figure 14 on the next page represents the time 
varying and time spread impulse response for two different values of cr3. 
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Figure 14. Time-Domain Representation of Doubly Spread Impulse Responses 

Based on the previous analysis such a channel is doubly spread and possesses 
both the time-variability and time-spread features. The effects of such a channel in the 
frequency domain are omitted because the analysis of the signal processing algorithms 
described in the next chapter is conducted only in the discrete time domain. 

C. LIMITATIONS AND CONCLUSIONS 

In this chapter, we have examined the time varying and spreading nature of the 
underwater medium. It was shown that the small-scale effects are significant when higher 
rate and coherent underwater communications are desirable. The static impulse response 
derived by a physics-based model like Bellhop was proven inadequate as an underwater 
channel simulation tool when the small-scale effects are important. Therefore, a new 
statistics-based simulation model was developed in this thesis in order to incorporate 
some of the small-scale effects of the real underwater medium. The new model is 
intended as a better simulation tool in order to test the receiver structures proposed in the 
chapters that follow. 


22 








The advantages of the proposed model are summarized below: 

• The impulse response is time variable. 

• Time spreading due to micropaths can easily be manipulated by 
varying the cr^ on the Gaussian. 

• Random phase shift was introduced in order to simulate micropath 
propagation 

• Spectral components are attenuated differently as a function of 
time-frequency. 

The proposed model has the following limitations: 

• It adds statistics on top of the statistics that were used in the first 
place to derive the Bellhop static impulse response. 

• The transmitted signal was an impulse and not a real waveform, 
which does not frilly develop the spectral distortion imposed by a 
real channel. 

• The degree of spread and micropath phase shift may be 
exaggerated compared to the real channel impulse response. 

In the next chapter we discuss the DSSS communication scheme and the proposed 
blind equalization algorithms to improve the receiver. 
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IV. ACOUSTIC SPREAD-SPECTRUM SIGNALING 


The previously proposed Seaweb signaling scheme in [1] is based on DSSS. It is 
DS-IQ-BPSK modulated at the transmitter and processed with a RAKE receiver. 
Although the transmitter proved to be adequate, the RAKE receiver did not prove to be as 
robust as anticipated. Therefore, a new receiver structure is needed for efficient 
equalization. This chapter briefly reviews the necessary theory for spread-spectrum and 
multi-rate signal processing and then analyzes the existing DSSS signaling scheme. 
Finally, this chapter examines and evaluates a new receiver structure based on multi-rate 
signal processing. The analysis emphasizes the signal processing aspect of the receiver 
and is confined to the discrete time domain. Several cases are evaluated for the blind 
equalizer based on whether receiver synchronization has been achieved. The actual 
implementation and integration in the existing DSSS scheme will be the task of future 
work. 

A. SPREAD-SPECTRUM AND MULTI-RATE SIGNAL PROCESSING 

I. Direct-Sequence Spread-Spectrum 

The choice of a DSSS transmission comes as no surprise, given its resistance to 
small-scale fading imposed on the received signal by the channel. This is achieved by 
spreading the signal energy over the widest available bandwidth using a Pseudo-Noise 
(PN) spreading sequence. This sequence possesses cross-correlation characteristics 
resembling white noise as much as possible. 

A common way to produce a PN sequence is by using an n-shift register, with the 
outputs of each stage combined by a network of exclusive-or (XOR) gates. The output of 
every XOR gate is fed back to the input of the shift register in an appropriate manner so 
that the PN generator will cycle through each of the possible states at a random order 
before any state is repeated. The resulting sequence is called a maximal-length sequence 
(or “m-sequence”). A serious disadvantage of m-sequences is that their cross-correlations 
can be quite large compared to the cross-correlation of true random binary sequences. 
However, a smaller subset of these m-sequences that possesses small cross-correlation 
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characteristics exists, called “preferred m-sequences.” A modulo-2 summation of a pair 
of these latter sequences results in a new pair of preferred m-sequences that possess the 
desired cross-correlation characteristics. Finally modulo-2 summing one of a pair of 
preferred m-sequences with ^-cyclically shifted versions of the other results in a Gold 
code. The interested reader can find more information on PN sequences and Gold codes 
in Ref. [12]. An example on the generation of a Gold code from two preferred m- 
sequences is illustrated in Figure 15 below: 



g'(D)=\+D^ + D^-¥D*-¥D^ 

Figure 15. An Example of a Gold-Code Generator [From Ref. 12.] 


A PN sequence has a fixed length of P chips. If we let be the pulse duration of 
the chipping sequence the PN sequence, is periodic with period P ■ . Furthermore, if 

is the signal pulse duration, the ratio of chips to signal is called the spreading gain k 
where 


k = P 



(4.1) 


Moreover, if we define the chipping rate as the inverse of and the symbol rate 
as the inverse of , we have the following relations: 


and 


Rr 




(4.2) 


(4.3) 


26 











Therefore, a chipping sequence of a chipping rate produces P = R^lR^ 


chips/symbol. If we combine the three previous relations we can write: 

k = P = ^ = -P. 


(4.4) 


The spreading is done at the baseband by multiplying the baseband data 
waveform with the PN sequence (also called the chipping sequence). In the receiver, the 
received signal is despread using a synchronized replica of the transmitter’s PN sequence. 
By using this technique, the channel induced interference is not despread along with the 
encoded data sequence because it possesses different correlation characteristics than the 
data sequence. Therefore the different multipath arrivals are rejected and the original 
signal is reconstructed as it “floats” over the noise. The effects of wideband interference 
on a DSSS signal are shown in Figure 16 below. 



Figure 16. Wideband Interference on the DSSS Signal [After Ref. 1.] 

A chipping sequence can be characterized as “short” when the entire sequence is 
transmitted within every data bit or “long” when only a portion of the sequence is 
transmitted within each bit. In the original Seaweb modem implementation, a long Gold 
code was used whereas in the proposed receiver structures a short code of length N 
proved to give better results. 
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2, Interpolation by an Integer Factor / (Upsampling by 1) 

In many practical applications, changing the sampling frequency of a given 
signal, either to increase it or decrease it, is desirable. The terms “interpolation” and 
“upsampling” are used interchangeably in signal processing to denote a desired increase 
in the sampling frequency , usually by an integer factor I. The analysis that follows is 
based on [13] and will generally retain the same symbols. 

If we let Fy represent the initial sampling frequency and Fy represent the increased 
sampling frequency, the upsampling operation can be expressed mathematically as: 

^ = I<^Fy=F,-I. (4.4) 

Fx 

Moreover, since the sampling interval is equal to the inverse of the sampling 
frequency or = 1/the operation can be rewritten in terms of the sampling interval as: 

= (4.5) 

where J),and Ty are the sampling intervals corresponding to the initial and to the 
increased sampling frequenciesand T].,respectively. In other words, this means that 
the sampling interval of the increased sequence should be / times smaller than the initial. 
This operation is usually done by adding 7-1 zeros between samples. Suppose the signal 
we desire to upsample is the discrete-time sequence x[n\. If we let the desired sequence at 
a higher sampling rate be y[n\, the upsampling process can be represented by the block 
diagram in Figure 17 below for the specific case of interpolating by three. 

T 


x[n\ y[n\ zeros 



T 


Figure 17. Upsampling by an Integer Factor I in the Time Domain 
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The upsampled sequence in the time domain retains all of the initial information. 
However, artifacts are generated by the additions of the zeros, which are evident in the 
frequency domain. Because the zeros were inserted in every/-1 sample the sequence 
y\n\ = y[m-I] has non-zero terms only in every m samples. The z-transform of the 
upsampled sequence is 

00 00 

7(z)= X y[m-I]z-’ = ^ x[mY”'^ = X{z‘) (4.6) 

m=—co'' ^ m— —CO 

x[m\ 

where v[w] is the initial sequence. If evaluated on the unit circle, the spectrum of the 
upsampled sequence F((yY)is 

Yio)y) = Xio)y-I) (4.7) 

where co^ is the radian frequency corresponding to the increased sampling frequency. 
Furthermore, the relation between co^ and cOy. (the initial radian frequency) is 


COy 


I 


(4.8) 


The effect of the upsampling operation in the Fourier domain is shown in Figure 18 on 
the next page for the particular case of upsampling by two: 
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Figure 18. Spectral Effect of Upsampling 

We observe that the initial spectrum is squeezed by / ( / = 2 in the case of Figure 
18). Therefore, if the spectrum of the initial sequence lay between-;r and ;r, the 
spectrum of the upsampled sequence is squeezed within the interval -Kj I and Kjl. As a 
consequence, image frequencies appear in the intervals[-;r,-;T//]and [;r//,;r] which are 

artifacts created by the process. Fortunately, the image frequencies can be easily 
eliminated with Low Pass Filtering after upsampling. 


3, Decimation by an Integer Factor D (Downsampling by D) 

Decimation is the inverse process of interpolation. In downsampling, we decrease 
the sampling frequency of a signal by an integer factor D. 


If we let be the initial sampling frequency and Fy be the decreased sampling 
frequency, the downsampling operation can be expressed analytically as: 


El. 

F, 


D 





(4.9) 
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or in terms of the sampling interval: 


1 


1 


1 1 


1 


Fy = -<:^>-=- =Ty - 

D Ty D D 


(4.10) 


where and Ty are the sampling intervals eorresponding to the initial and the deereased 
sampling Ifequeneies F^ and Fy, respeetively. The most usual representation is shown in 
Figure 19 below for the speeifie ease of deeimating by three: 



Figure 19. Downsampling by an Integer Faetor D in the Time Domain 

In the time domain deeimation ean be thought of as retaining one for every D-\ 
samples and information is irretrievably lost during this proeess. 

As expeeted, downsampling ereates aliasing in the Ifequeney domain. In faet Ref. 
[13] shows that the downsampled sequenee has a eorresponding speetmm Y{cOy)3i's, 
shown in the following equation: 


Dto [d D, 


(4.11) 


where cOy is the radian frequeney eorresponding to the deereased sampling ffequeney. 
Moreover, the relation between the inereased radian Ifequeney cOy and the initial radian 
frequeney co^ is 


cOy=D-cOy^. (4.12) 

The effeet of the downsampling operation in the Ifequeney domain is shown in Figure 20 
below for the partieular ease of downsampling by two: 
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Figure 20. Spectral Effect of Downsampling 

We observe that the initial spectrum is stretched by D. Therefore, if the spectrum 
of the initial sequence lay within the interval -;r/D and ;r/D, the spectrum of the 
downsampled sequence is stretched within the interval-;r to n. The immediate 
consequence is that frequencies below -;r/Dor above ;r/D of the initial sequence are 
folded back and aliased. In Figure 20, aliases were eliminated by low-pass filtering the 
original sequence in the interval [-;r/2,;r/2] before downsampling. 


4, Upsampling and Spreading 

As seen previously, the spreading of a signal is done at the baseband by 
multiplying the baseband data waveform with the PN sequence [12]. Since all of the 
analyses that follow are done in the discrete time domain and the signals involved are 
impulses, the effect of spreading the data sequence is equivalent to multiplying each bit 
of the data sequence with the PN sequence. From a signal processing point of view, this 
is equivalent to upsampling by P and then convolving with the spreading code. The scope 
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of the following Figure 21 is to demonstrate the equivalence between the two claims for a 
single bit. 
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Figure 21. Upsampling and Spreading a Sequence 


As we can observe, the result from both methods is identical. Following the 
upsampling with zeros and convolution approach, each bit is appended with the PN 
sequence. Equivalently following the multiplication approach, each bit appears to be 
spread throughout the PN sequence. 


B. PREVIOUSLY PROPOSED DSSS SIGNALING SCHEME 

In the previously proposed DSSS signaling scheme, a DSSS technique is used to 
chip and to spread a data bit sequence, which has been modulated by a Differentially- 
encoded Binary Phase-Shift Keying (DBPSK) carrier. The process is depicted in Figure 
22 on the next page. 
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Figure 22. Previously Proposed DSSS Scheme [After Ref. 1.] 

The transmitter scheme was implemented as part of a previous thesis [1]. The 
main problem in the previous implementation was the RAKE receiver. Although this type 
of receiver is used in cellular communications, it did not perform as anticipated 
underwater. 

The RAKE receiver is based on the principle that multipath arrivals separated by 
time intervals greater than the coherence time of the channel T^can be resolved and 

summed. This technique is also referred to as the “integrate-and-dump” method and the 
interested reader can find more details in [11]. There are two major problems associated 
with this technique. First, there is the need to estimate the amplitudes of the equalizer 
continuously over a time-varying channel which incurs a great computational expense. 
Moreover, in the case that the channel is fast fading, there is no guarantee that the 
equalizer will converge before the channel changes again. Secondly, when in any 
processing leg of the structure there is no signal present, the input is due to the noise 
only. However, the latter is resolved by selecting an appropriate threshold. The complete 
structure of the RAKE receiver is shown in Figure 23 on the next page. 
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Vt 

Figure 23. Block Diagram of a RAKE Receiver for DS-IQ-DBPSK [From Ref. 1.] 

The present thesis follows a different approach in the equalization process, in 
which equalization is obtained by subspace-decomposition. 

C. PROPOSED DSSS SCHEME 

1. Simplified Transmitter Structure 

Two blind equalization algorithms, based on multi-rate signal processing and 
subspace-decomposition for DS-CDMA, are evaluated in this thesis for various 
interference and synchronization assumptions. For the sake of simplicity, the transmitted 
signal is assumed to be a bit sequence, which is then chipped and transmitted through the 
channel. All of the analysis that follows is done in the discrete time domain. Channel 
encoding and modulation are not implemented in this thesis, but they are assumed to 
improve the signaling scheme. Incorporating the simplified transmitter and receiver 
structures used in this thesis in the actual Seaweb scheme will be the task of future work. 
The simplified transmitter used throughout the remainder of the thesis is as shown in the 
following block diagram in Figure 24 on the next page. 
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Figure 24. Simplified Multi-Rate DSSS Transmitter 
The data sequence a[n\ is first upsampled by P = T^jT(. and then chipped by the 
chipping sequence c[n]. The chipping sequence is a short PN code of length P. 
Essentially, this means filtering the upsampled data sequence a[n] by the P-th order 
filter: 

C(z) = c[0] + c[l]z-‘ +... + c[P - l]z-<'’-‘^. (4.13) 

The state space equivalent transmitter structure is represented in Figure 25 below: 


C(z) 



Figure 25. State-Space Representation of CDMA Transmitter 


Subsequently, the chipped signal is convolved with the channel impulse response 
g[n\. The channel is assumed to be of length M < P but, as will be shown in the next 
chapter, the multi-rate technique was proven robust even when the whole impulse 
response of the underwater channel was used. The corruption by the channel is equivalent 
to filtering by the M-th order linear filter: 

G{z) = g[0] + g[l]z-‘ +... + g[M - l]z-<^-'>. 
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(4.14) 






















Finally, noise is added to account for random interference. As we will see later in 
the chapter, the noise will be both Gaussian white noise and colored noise, in order to 
evaluate the receiver’s behavior in a simulated adverse environment before using the 
time-variable impulse response. 

2, Blind Equalizer Structure and Matched Filtering 

The receiver structures that follow are considered blind equalization methods. By 
the term “blind equalization,” we mean that the receiver has no a priori knowledge of 
what the transmitted sequence nor the corrupting channel is but attempts to undo the 
channel effect, as well as adjusts the fdter coefficients based on statistical methods 
applied to the received signal. The receiver structure is based on a DS-CDMA algorithm 
from Ref. [14] and is adopted for a single user. The matched filter analysis has been 
adapted from Ref. [15]. Such algorithms have been described for cellular 
communications and are simulated in the current thesis for the underwater medium. The 
blind equalization algorithms presented in this thesis are intended for off-line use at the 
present stage since the entire received sequence must be available in order for the 
algorithms to work. Moreover, the synchronization between the receiver and the 
transmitter is assumed to be achieved using the current Seaweb synchronization methods. 

a. Additive Gaussian White Noise (AWGN) without Multipath 

In this case, the g[n] portion of the channel is omitted from the transmitter 
in Figure 24. The data sequence is upsampled to the chip rate by P and then multiplied by 
c[n], which has the same effect as appending a code to each bit. The resulting sequence is 
transmitted and Gaussian noise is added. Therefore, the received sequence y[n] is simply 
the chipped sequence with added Gaussian white noise and the resulting vector length 
does not change. 

The receiver structure for the white noise case is shown in Figure 26 on 
the next page. The received sequence y[n] is subdivided into P vectors of length equal to 
the original data sequence. This is achieved by downsampling the original sequence by P 
to obtain the first vector, then forward-delaying by one sample and downsampling by P 
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to obtain the next vector. The process is repeated P times in total. The optimum filter is 
derived from the resulting matrix . 


An] 



a\n\ 


Figure 26. Multi-Rate CDMA Receiver for AWGN 

To gain more insight for the algorithm, we combine the transmitter and the 
receiver in the manner shown in Figure 27 below: 



Figure 27. Combined Transmitter-Receiver for AWGN 

The receiver structure seems non-causal as it forward-delays and 
downsamples the received sequence. However, this is not the case. Suppose we select the 
k-th branch from the transmitter part and the ^-th branch from the receiver (where k, q= 
0, 1, ..., P -1) as shown in Figure 28 on the next page: 
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yk.M] 


a[n] 



> 


Figure 28. Multi-Rate Representation of a Branch 

If we combine the two delay terms into one and observe that there is no 
cross-talk the branch is simplified to the following Figure 29: 



Figure 29. Simplified Multi-Rate Representation of a Branch 

f 0, k ^ q 

Since there is no cross-talk, yt,\n\ = { and the 

\c{k]a{n\, k = q 

q-k exponent reduces to zero when k^q. Therefore, the branch is causal and 
Ft <,[”] reduces to yiXn\ = c{k~\a\n\ + w^\n\. Given that, for every k = q only one branch 
is non-zero, we can generalize for the whole system: 

y[n] = c • a[n] -I- w[n] (4.15) 



y[nP] 


■ c[0] ■ 

where the vectors y[n] and c are defined as: y{n\ = 

y\nP + P - \] 

and c = 

c[P-\]_ 


Using the previous equation and defining a{n\ to be the estimate of a{n\ 
(the transmitted signal) then 

«[«] = f yin] = f{c■a{n] + ^^n]) (4.16) 

where is the desired filter (estimate of the channel), as shown in Figure 26, which is 
yet to be computed. The best estimate for / is the one that maximizes the SNR in 
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analogy with the matched filter. The matched filter is described by the following 
equation: 


{\fc-a[n]f] 


(4.17) 


If we assume a\n\ to be stationary, and a{n\ and iy[n] independent, 

^ {l a[n\ f } reduces to . Furthermore, since £■ {| fv^n\ f J = f R^f_ = f <7J ■ f = 
cr^f^f, the SNR can be rewritten: 


SNR = 


crjff 


(4.18) 


The SNR is maximized when f = A-c with A any scalar constant. 
Therefore, the desired filter in the white noise case is 

f = A-c. (4.19) 


b. Multipath with A WGN 

The transmitter block diagram for the multipath and additive white-noise 
case is as shown in Figure 24. The data sequence is upsampled to the chip rate by P, 
spread by a short PN sequence and then convolved (corrupted) with the channel impulse 
response g\n \. Finally, white noise is added to compensate for random interference. The 
receiver structure is slightly different than the previous case. The received sequence is 
first despread using a synchronized replica at the receiver and then forward delayed and 
downsampled .The previous analysis on the causality of this process is still valid and will 
not be repeated. The receiver block diagram is shown in Figure 30 on the next page. 
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Figure 30. Receiver Structure for Multipath with AWGN 
In Figure 30 above the despreading sequence shown as c[-n\ is such that 
c[n\*c[-n\ = 5[n\. The combination of transmitter-receiver is shown in Figure 31 below: 



Figure 31. Combined Transmitter-Receiver for Multipath with AWGN 


Following the previously laid methodology for the white noise case the 
matrix y[n\ is composed by P versions of the original data sequence, which has been 

corrupted by the channel. This can mathematically be expressed as: 

An] = g-a[n] + )An] (4.20) 

where g is the unknown channel corruption and w[n] the additive white noise. Since the 
equalization method is blind, the signal sequence is also unknown. However, we can 
make an estimate of a{n\ called a{n\ = y{n\ = (g •«[«] +provided we can 
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estimate /^(the desired fdter). Following the same methodology as previously we will 
seleet the that maximi z es the SNR in analogy with the matehed fdter. The matehed 

filter, however, presupposes knowledge of the transmitted signal, whieh in our case is 
unknown. To overcome that, we use a different metric of SNR called SNR+i [16] in order 
to make use of the received sequence y[n]. The new matched fdter is described by the 
following equation: 


SNK, = ■ 


E\\fg-a[n-\\ 


+ 1 = 


E{\fy[n]/[n]f\} 


fRJ_ 


{l/VMf} E^f^^n\v/{n\f^^ fR^f_ 


(4.21) 


Since f = f crj ■ f = cr^f f_ the SNR^^ can be rewritten: 


SNK, 


fRJ_ 


(4.22) 


In order to maximize the SNR^^wq set the denominator to one and 
maximize the numerator. Therefore, we choose / so that it maximizes the quantity: 

2. = rRJ_ + X{\-CTyf_) (4.23) 


where T is a real Lagrange multiplier. The maximum for R^ can be obtained by taking 
the gradient with respect to / and setting it to zero: 

V, L = RJ_ - XaJ « R^f = Aajf « R^f = A'f , (4.24) 

which is an ordinary eigenvalue problem and / is the eigenvector that corresponds to the 
maximum eigenvalue the autocorrelation matrix Ryy = E^y[n\y^[n^. Applying 

the matched fdter to the received signal yields: 

a[n] = {fYy[n]. (4.25) 
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c. Multipath with Additive Colored Noise 

In the case that the added noise is colored, the receiver would still have the 
same structure as previously (Figures 30 and 31) but the derivation of the optimum filter 
/ would have to account for the fact that the noise is not white but colored. In this case 

the autocorrelation matrix of the noise is not diagonal and the off-diagonal terms will 
possess non-zero values. If we let the new autocorrelation matrix for colored 

noise, the eigenvalue decomposition of would yield: 

R =E A (4.26) 

WW WWW V / 

where E^ is the eigenvector matrix and is the eigenvalue matrix of R^^. Using the 
Mahalanobis transformation: 

y [n] = R~\^-y[n] = (E^AjEj)-y[n], (4.27) 

the autocorrelation matrix A of the colored noise is transformed into a new coordinate 
system in which the noise can be considered white. Applying the coordinate 
transformation (4.23) on the original autocorrelation matrix R^ yields: 

R'yy=R^\^RyyR^^^- (4-28) 


Then applying Equation (4.21) to the new coordinate system, we can 

write: 

R'„f' = ^f', (4.29) 

which essentially means that the problem reduces to a regular eigenvector problem in the 
new coordinate system and /' is the eigenvector that corresponds to the maximum 

eigenvalue of the autocorrelation matrix R'^ on the new coordinate system. Applying 
the matched fdter in the original coordinate system is done by: 

a[n] = (/') y {n\ = (/') R^ ^ y[n\ = f y[n] 

where the optimum filter in the old coordinate system is 
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(4.30) 



(4.31) 


f={f'^ R "2 

J \^ / WW 


d. Loss of Receiver Synchronization 

The foundation for the formerly developed receiver algorithm is that in 
order to remove the channel effect adequately, synchronization between the transmitter 
and receiver is imperative. Due to synchronization there is only one dominant 
eigenvector of the autocorrelation matrix which belongs to the largest eigenvalue. 

However, synchronization might not always be possible in the adverse underwater 
medium. In the present analysis, we will examine a possible remedy for the loss of 
synchronization. 

The basic receiver structure is the same as in Figure 30. As shown 
previously, upsampling by P (at the symbol rate) in the transmitter extended the bit 
through the whole length of the code. Therefore, downsampling and forward delaying at 
the receiver results in the matrix y[n], which comprises the P versions of the original data 

sequence a[n], corrupted by the channel. The current case is diversified from the previous 
because the loss of synchronization in the PN sequence causes the receiver to despread 
the received sequence incorrectly as shown in the following Figure 32: 


Receiver Synchronized with Transmitter 


Transmitter 

Receiver 

Transmitter 

Receiver 


Receiver Not Synchronized with Transmitter 


Despread 
in error 



Obscured by the Unobscured Obscured by the 

previous bit portion of bit next bit 


Figure 32. Loss of Receiver Synchronization 
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Nevertheless, it is shown in the following analysis that the loss of 
synchronization is an anomaly that appears in the eigenvalues of the autocorrelation 
matrix and spans a limited subspace. Let us assume that each portion of y[n\ contains 

two symbols at maximum, namely the current and a portion of the previous. This is 
equivalent in writing to: 

An\ = g^a[n'\ + g^a[n-\'\ (4.32) 


where g^ and gj depend on both the channel and the autocorrelation of the code. In 
matrix form, (4.32) can be expressed as: 


An] = \ga,gy 


a[n] 

a[n-l] 


(4.33) 


From Equation (4.33) we can conclude that y\n\ spans a two-dimensional 
subspace defined by the vectors and gj. In order to estimate the new subspace, we 
will examine the eigenvalues and eigenvectors of the autocorrelation matrix Figure 
33 below illustrates an example of what the eigenvalues of i?^look like when there is no 
synchronization: 



Figure 33. Plot of the Normalized Eigenvalues When No Syncronization Is Achieved 
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The example uses a particular case with upsampling P = 50 and the 
MATLAB sorting of the eigenvalues in which the largest eigenvalue is placed last. We 
observe that only two eigenvalues are significant whereas the others are very close to 
zero and can be neglected. 

In the absence of noise, this means that y[ri\ can be analyzed by these two 
nonzero eigenvalues and corresponding eigenvectors y and only: 


y[n] = ye^[n] + y2e^[n]. 


Since the eigenvectors are orthonormal: 


e\n] = y y[n] and e^in] = y^y[n\ ^ e[n] 



Substituting (4.29) into (4.31) yields: 

a\n\ 
a[n -1] 



(4.34) 


(4.35) 


(4.36) 


Assuming and gj are linearly independent and if we let 


T 

f 

T 




T 

1 1 

T 

— 

T 


nil 

\ 


y 


can write: 


a[n] 


T 

a\n-\ 


T 

^2 


■e[n\. 


(4.37) 


Since Si ^^e linearly independent, it follows that Wj^and are 

also linearly independent. Consequently, since a[n] = ■ e[n] it is true that 

a[n -1] = ■ e[n] a{n\ = ■ e[n + 1] .In vector form the system of equations can be 

expressed as: 


1 T T\ 

e[n\ 

m, ,-^2 • 

e[n + 1] 


0 . 


(4.38) 
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In order to calculate m 


-nil 


, we start by calculating the covariance 


matrix = 


e[n] 
e[n + 1] 


e^[n],e^[n + l], then estimate the eigenvalues and eigenvectors 


and select the eigenvector corresponding to the minimum eigenvalue to be w. The 
estimated sequence a[n\ is derived from (4.37), once the estimates for m are known. 


D. SUMMARY 

In the current chapter, we first examined the theory relevant to DSSS and multi¬ 
rate signal processing in order to better appreciate the previously proposed DSSS 
communication scheme and the proposed receiver changes. Subsequently, the problems 
inherent in the previously proposed DSSS communication scheme were shortly reviewed. 
Although the transmitter scheme was proved to be sufficient in the original Seaweb 
communication scheme, the adaptive implementation of the RAKE receiver was 
problematic. Finally, the current chapter proposes a RAKE receiver replacement. Two 
blind equalization algorithms were developed in order to circumvent the adaptive tap 
estimation, based on whether the synchronization between the receiver and the 
transmitter was achieved. Both algorithms are based on the fact that the channel 
interference will always span a small subspace in the eigenvalues of the autocorrelation 
matrix of the received signal. If the subspace is very large, that means that the time 
coherence of the channel is too small, which leaves no margin for coherent 
communications of any kind. Both algorithms were analyzed in the discrete time domain 
for the various combinations of multipath, white and colored noise. The simulation of 
both algorithms in the various combinations of noise and multipath was done in 
MATLAB and the results are presented in the next chapter. 
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V. IMPLEMENTATION AND PERFORMANCE 


This chapter discusses the MATLAB implementation of the proposed Seaweb 
receiver design and the performance results in simulated ocean channels. The chapter 
starts with the necessary reference to Ef^jN^ and Bit-Error Rate (BER) and then presents 

the performance of the receiver in terms of BER vs.E^INq for the various combinations 

of multipath and noise conditions. The block diagrams of the transmitter and the receiver 
structures were presented in Figures 24 through 31 of Chapter IV for the various 
combinations of multipath and noise. Each case is discussed here separately with 
particular emphasis on how varying specific design parameters (e.g., bit rate, packet size 
etc.) affects the (BER) of the receiver. Under the current configuration Seaweb operates 
with 72-bit utility packets and an available bandwidth of 5 kHz. The utility packet 
information bit-rate when transmitted with DSSS is proposed to be 40 to 100 bps. Finally, 
the performance of the receiver is measured with a simulated doubly spread underwater 
acoustic model and additive colored noise case, which is a representative case of the 
underwater channel. 


A. PERFORMANCE ANALYSIS OF COMMUNICATION SYSTEMS 
I. Signal Power (S) 

The power of the transmitted signal S is evaluated using the expression: 


S = ^ - 

NE 


(5.1) 


where x[n\ are the amplitudes of the DSSS signal samples and N is the total number of 
chips. If we normalize so that |x[i?]| = 1, (5.1) can be rewritten as: 


5 


A^-l 

III 


n=0 


N 



(5.2) 
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2, Signal-to-Noise Ratio (SNR) 

The SNR of the received signal was defined in Equation (2.2) as the ratio of the 
time average signal power divided by the time average noise power. Practically the SNR 
is the signal power excess over the channel noise and is a natural figure of merit used to 
describe the performance of communications. In digital communications, it is more 
common to use a scaled version of the SNR called E'j jN^ where is the energy-per-bit 
and is the noise power spectral density. The analytical expression for the relation 
between E^/N^ and SNR is the following: 

E W 

-^ = SNR - (5.3) 

No R, 

where W is the available bandwidth of 5 kHz and R^the bit rate. It is mentioned in [I] 
that when the DSSS signal SNR < -6 dB the signal cannot be audibly distinguished from 
the background noise, however it is loud enough to be detectable by the receiver. For a 
40-bps signal the corresponding E^jN^ threshold is approximately22 dB. In order to 
enable a direct comparison with [1], the 22 dB threshold was maintained in this thesis. 
The range of the desired E^^jN^ used for the simulation was from 0 through 25 dB. 


3. Bit-Error Rate (BER) 

BER is the metric that quantifies the reliability of a communication system as the 
ratio of the erratic bits to the total bits sent: 


BER = 


Errors 

Total Number of Bits 


(5.4) 


It depends on the ratio E'^/A^o their interrelation is important when a small 
E^/Nq must be maintained while the BER must be kept to a minimum (below 10“^ for 
Seaweb). Therefore, the performance throughout the simulation was expressed in terms 
ofBER vs.E^/A^g plots. 
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B. PERFORMANCE RESULTS 


The underwater aeoustie ehannel was modelled as deseribed in Chapter III in 
order to ineorporate the small-seale effeets. However, the performanee of the reeeiver 
was measured for all of the eases examined in Chapter IV before being measured for the 
doubly spread ease. The simulation involved sending a thousand packets for each of the 
values in order to create sufficient statistics. The BER plots therefore represent 

Monte Carlo realizations of the performance in the various conditions. 

1. AWGN Only 

The case of the AWGN channel was chosen in order to test the equalization 
algorithm and obtain an estimate of the receiver performance in the absence of multipath. 
In the case of AWGN channel, white noise is added to the transmitted signal and the 
received signal is of the formy[n] = c[n\ * a[n\ + rv[n]. The receiver configuration used is 
that described by Equations (4.15) through (4.19) and represented in Figure 26. 

Figure 34 below demonstrates the receiver performance in AWGN when the 
packet size is varied from the required 72 to 1242 bits while the bit-rate is fixed at 40 bps. 
We observe that increasing the packet size does not degrade the receiver performance. 



Figure 34. Performance in AWGN Only When Varying the Packet Size and Bit-Rate is 

Fixed at 40 bps 
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2, Static Channel Only 

In this case, the multipath effect is isolated and the reeeiver performanee is 
evaluated for stability over 100 trials. The unfiltered statie impulse response, as shown in 
Figure 12, was used for this ease. However, sinee there is no noise available, there ean be 
no BER ys.Ei^INq evaluation. The error rate over the number of trials was used as a 

performanee metrie in this ease. Sinee the impulse response is statie, the reeeived signal 
is of the formy[n\ = g[n\*{a[n\*c[n\). The reeeiver eonfiguration evaluated is that of 
Figure 30, and the reeeiver sueeessfully equalized the multipath effeet in all of the trials. 

3. Static Channel with AWGN 

This is the natural eombination of the two previous eases in whieh the signal is 
subjeeted to the eombined effeet of multipath and AWGN. The reeeived signal has the 
form y[n\ = g[n \* a\n\ + >v[«] and the BER vs. is evaluated when different design 

parameters are varied. The reeeiver strueture evaluated is that of Figure 30 deseribed by 
Equations (4.20) through (4.25), and the reeeiver performanee is demonstrated in Figures 
35 through 40 on the following pages. 

Figure 35 demonstrates the reeeiver performance when the paeket size is varied 
from the required 72 to 1242 bits while the bit-rate is fixed at 40 bps. As expeeted from 
the previous eases examined so far, inereasing the paeket size does not degrade the 
reeeiver performanee. 

Figures 36 and 37 demonstrate the reeeiver performanee for 72 and 144-bit 
paekets when the bit rate is varied in the range 40 to 200 bps. It is shown that the highest 
aehievable bit-rate for the settings used is 100 bps. 

Finally, Figures 38 through 40 show the effeet of the ehannel length on paekets of 
sizes 72 through 1242 bits at a fixed bit-rate of 40 bps. This ease is interesting beeause up 
to now we eluded the faet that the algorithm is supposed to work up to a ehannel length 
of P. This length ean be eonsidered as the eoherenee time of the ehannel. However, it is 
demonstrated in the aforementioned figures that the algorithm proved robust for mueh 
longer impulse response lengths. The ehannel lengths are stated in samples. 
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Figure 35. Performance in Multipath with AWGN When Varying the Packet Size and Bit- 

Rate is Fixed at 40 bps 



Figure 36. Performance in Multipath with AWGN When Varying the Bit-Rate and Packet 

Size is Fixed at 72 bits 
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Figure 37. Performance in Multipath with AWGN When Varying the Bit-Rate and Packet 

Size is Fixed at 144 bits 



Figure 38. Performance in Multipath with AWGN When Varying the Channel Length (in 
samples) for Bit-Rate 40 bps and Packet Size 72 bits 
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Figure 39. Performance in Multipath with AWGN When Varying the Channel Length for 

Bit-Rate 40 bps and Packet Size 144 bits 



Figure 40. Performance in Multipath with AWGN When Varying the Channel Length for 

Bit-Rate 40 bps and Packet Size 1242 bits 
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4, Colored Noise Only 

The colored noise was produced by filtering white noise by a 7* order FIR filter 
as follows: 


H{z) 


_ 1 _ 

z® - 0.9z® + 0.4z"* - 0.2z^ + 0. Iz^ - 0. Iz 


(5.5) 


This manipulation changes the form of the covariance matrix of white noise 

from , to non-diagonal by introducing off-diagonal terms so that the receiver could 

be tested for colored noise. In Figure 41 below white noise is shown to equally affect the 
frequency spectrum whereas colored noise is frequency-dependent. 



0 2 4 6 8 10 12 14 16 18 20 

Frequency (kHz) 



Frequency (kHz) 

Figure 41. Comparison of White and Colored Noise Effect in the Frequency Domain 

Figure 42 on the next page demonstrates the receiver performance when the 
packet size is varied from the required 72 to 1242 bits while the bit-rate is fixed at 40 bps. 
As expected from the previous cases examined so far, increasing the packet size does not 
degrade the receiver performance. 
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Figure 42. Performance in Colored Noise Only When Varying the Packet Size and Bit- 

Rate is Fixed at 40 bps 

5, Static Channel with Additive Colored Noise 

This case combines the effect of multipath and colored noise. The received signal 
has the form y[n\ = g{ri\ * {a[n\ * c[n\) + n{n\ and the BER vs.E^IN q is evaluated when 
different design parameters are varied. The receiver structure evaluated is that of Figure 
30, adapted for colored noise and the receiver performance is shown in Figures 43 
through 48 on the following pages . 

Figure 43 demonstrates the receiver performance when the packet size is varied 
from the required 72 to 1242 bits while the bit-rate is fixed at 40 bps. The simulation 
results are consistent with the previous cases examined so far and increasing the packet 
size does not degrade the receiver performance. 

Figures 44 and 45 demonstrate the receiver performance for 72 and 144-bit 
packets when the bit rate is varied in the range 40 to 200 bps. It is shown that the highest 
achievable bit-rate for the settings used is 100 bps. 

Finally, Figures 46 through 48 show the effect of the channel length on packets of 
sizes 72 through 1242 bits at a fixed bit-rate of 40 bps. As previously, the algorithm 
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proved robust for impulse response lengths longer than P. The channel lengths are stated 
in samples. 



Figure 43. Performance in Multipath with Additive Colored Noise When Varying the 
Packet Size and Bit-Rate is Fixed at 40 bps 



Figure 44. Performance in Multipath with Additive Colored Noise When Varying the Bit- 

Rate and Packet Size is Fixed at 72 bits 
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Figure 45. Performance in Multipath with Additive Colored Noise When Varying the Bit- 

Rate and Packet Size is Fixed at 144 bits 



Figure 46. Performance in Multipath with Additive Colored Noise When Varying the 
Channel Length (in samples) for Bit-Rate 40 bps and Packet Size 72 bits 
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Figure 47. Performance in Multipath with Additive Colored Noise When Varying the 
Channel Length (in samples) for Bit-Rate 40 bps and Packet Size 144 bits 



Figure 48. Performance in Multipath with Additive Colored Noise When Varying the 
Channel Length (in samples) for Bit-Rate 40 bps and Packet Size 1242 bits 
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6, Doubly Spread Channel Only 

In this case the effeet of multipath in a doubly spread channel is isolated from 
noise and the performanee of the reeeiver eonfiguration in Figure 30 is evaluated for 
stability. The impulse response shown forcr^ = 8 in Figure 14 was used in this case. Sinee 
there is no noise available, there can be no BER vs.E’^/A^o evaluation. The percentage 

error rate over the channel length was used as a performance metrie in this case, given the 
faet that the equalizer is designed to equalize an interferenee of length M <P . Sinee the 
impulse response is time-variable, the ehannel is considered statie only within the 
eoherenee time of the channel. The reeeived signal can only be obtained by eonvolving 
the transmitted signal with the channel, if the signal duration is eomparable to the 
eoherenee time of the channel. The performanee of the receiver in the doubly spread case 
is shown in Figure 49 below for the speeifie case where P = 50. We observe that, 
although the algorithm outperformed equalization up to a channel length of 3'P = 150, 
the error rate inereases abruptly to 50% for longer ehannel lengths. 
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Figure 49. Error Rate vs.Channel Length (in samples) for Doubly Spread Channel 
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7. Doubly Spread Channel with Additive Colored Noise 

This case was chosen to test the performance of the receiver shown in Figure 30 
for the adverse case of a doubly spread channel with additive colored noise. Since the 
channel impulse response is time-variable, the received sequence can only by derived 
from the transmitted sequence by convolution with the channel impulse response if the 
signal duration is comparable to the coherence time of the channel. The colored noise is 
additive to enable the BER vs. evaluation. The receiver performance is 

demonstrated in Figures 50 through 55 on the following pages. 

Figure 50 demonstrates the receiver performance when the packet size is varied 
from the required 72 to 1242 bits while the bit-rate is fixed at 40 bps. The simulation 
results are compatible with the previous cases and increasing the packet size does not 
degrade the receiver performance; however, more power is needed to achieve the desired 
BER. 

Figures 51 and 52 demonstrate the receiver performance for 72 and 144-bit 
packets when the bit rate is varied in the range 40 to 100 bps. It is shown that the highest 
achievable bit-rate for the settings used is 40 bps. 

Finally, Figures 53 through 55 show the effect of the channel length on packets of 
sizes 72 through 1242 bits at a fixed bit-rate of 40 bps. It was shown previously that in 
the absence of noise the maximum channel length that can be equalized using this method 
is 3 - P = 150 (samples). Therefore the channel lengths simulated here were in the range 
50 to 150 (in samples). Consistent with the previous findings, the algorithm proved 
robust for impulse response lengths longer than P. 
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Figure 50. Performance in Doubly Spread Channel with Additive Colored Noise When 
Varying the Packet Size and Bit-Rate is Fixed at 40 bps 



Figure 51. Performance in Doubly Spread Channel with Additive Colored Noise When 
Varying the Bit-Rate and Packet Size is Fixed at 72 bits 
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Figure 52. Performance in Doubly Spread Channel with Additive Colored Noise When 
Varying the Bit-Rate and Packet Size is Fixed at 144 bits 



Figure 53. Performance in Doubly Spread Channel with Additive Colored Noise When 
Varying the Channel Length (in samples) for Bit-Rate 40 bps and Packet Size 72 bits 
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Figure 54. Performance in Doubly Spread Channel with Additive Colored Noise When 
Varying the Channel Length (in samples) for Bit-Rate 40 bps and Packet Size 144 bits 



Figure 55. Performance in Doubly Spread Channel with Additive Colored Noise When 
Varying the Channel Length (in samples) for Bit-Rate 40 bps and Packet Size 1242 bits 
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8. Loss of Receiver Synchronization 

In Seaweb, the receiver synchronization is achieved by appending three short 
pulses for acquisition and an FM chirp in the beginning of each data packet. When the 
PN sequence in the transmitter and the receiver are synchronized, the received DSSS 
sequence is despread by simple multiplication with the PN sequence in the receiver. 
However, due to the volatility of the underwater medium it is not unlikely that the 
acquisition frame may be lost, and the whole packet has to be retransmitted because the 
packet cannot be despread. It was shown earlier in Chapter IV that in the case where 
synchronization is lost, the resulting interference will appear in the eigenvalues of the 
covariance matrix of the received signal. Therefore, the recovery of the packet without 
retransmission is possible by adapting the receiver structure of Figure 30 according to 
Equations (4.32) through (4.39). 

Loss of synchronization is represented in Figure 32. In order to simulate loss of 
synchronization, a random delay of duration up to half a bit was injected to the received 
signal so that each bit contains a portion of the previous and the next bit. The receiver 
performance was only tested for the doubly spread channel and additive colored noise 
case. The algorithm switches to loss-of-synchronization mode when the second largest 
eigenvalue of the covariance matrix of the received signal is greater than two tenths of 
the largest eigenvalue. The receiver performance when the packet size is varied from 72 
to 1242 bits, while the bit-rate remains constant at 40 bps, is shown in Figure 56 on the 
next page. It is shown that the algorithm can achieve a BER of 10"^ or less for the 
specified boundary of 22 dB depending on the packet size used. The rest of the results 
closely match the results of the previous case when synchronization was achieved and 
will not be repeated. 
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Figure 56. Performance in Doubly Spread Channel with Additive Colored Noise When 
Synchronization is Lost as a Function of Packet Size and Fixed Bit-Rate at 40 bps 

C. SUMMARY 

In the current chapter, we first reviewed the relevant theory for communication- 
systems performance analysis and then implemented the receiver structures that were 
presented in chapter IV using MATLAB. The receiver performance was presented either 
as a function of BER ys. or channel error rate. Each case was simulated by 

varying several parameters (i.e., bit-rate, packet length and channel length) independently 
and measuring their effect on the receiver performance. Monte Carlo simulation was used 
for all simulated cases and a thousand packets were sent for each / Aq value in order to 

gather sufficient statistics. All of the algorithms proved to work adequately even for the 
adverse case of the doubly spread channel with additive colored noise. In most cases, the 
BER reached 10“^ below the E^jN^ =22 dB required for the system to possess LPD 
characteristics and under no conditions was there a need for retransmission due to 
excessive errors. 

In the next and final chapter, we summarize the findings and present areas for 
future work. 
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VI. CONCLUSION 


The goal of this thesis was to study a RAKE receiver replacement for use in 
Seaweb. The first step in the current research was to simulate a more realistic channel 
with time variability and time spread characteristic of undersea acoustic propagation. 
Therefore, a statistics-based model was implemented first in order to incorporate the 
small-scale effects inherent in the ocean. Subsequently, two blind equalization algorithms 
based on subspace-decomposition were implemented for use when receiver 
synchronization could be achieved or not. The performance of the receiver structures was 
first evaluated using a static impulse response for the various combinations of multipath 
and colored or white noise, and then for a simulated doubly spread channel with additive 
colored noise. A Monte Carlo simulation was used to measure performance in all cases, 
and the algorithms proved to work adequately in all of the cases simulated. 

A. FINDINGS 

The receiver structures simulated were shown to equalize the channel effect 
adequately, in most cases below the E'j/Ao=22dB limit required for LPD 

communications. A BER of 10~Ms generally possible even if that must be done at the 
expense of power. In the event that the 22-dB is a stringent boundary, a BER 

exceeding the 10"^ required for Seaweb was possible and errors could be mitigated 
through forward-error correction or automatic repeat request. 

The use of a 144-bit packet combined with a higher bit-rate of 100 bps was 
proven to be possible, except when the channel is very severe. As a comparison with the 
existing scheme, the packet size and bit-rate are twice or more than that required. In an 
extremely adverse environment, like the one simulated in the doubly spread channel with 
additive colored noise case, the higher bit rate of 100 bps was not possible. In any case, 
communications were achievable for the larger packet of 144 bits even if the bit-rate had 
to slow from 100 bps to 40 bps. 
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The modelled doubly spread channel allowed the testing of the algorithm in an 
environment in which small-scale phenomena like the coherence time of the channel are 
important. Given that the packet duration used was sufficiently small compared to the 
channel (except for the 1242-bit packet), we were able to test the algorithm in a time- 
varying and time-spread environment, which resembles the real ocean. The algorithm 
was successful up to a length three times longer than the specifications of M < P in all 
combinations of bit-rate and packet size. 

The loss of receiver synchronization algorithm was tested for a short PN sequence 
and for a shift up to half a bit and was proven to work with the settings shown. A BER of 
10"^ or less for 22 dB was proven possible in a simulated doubly spread 

channel with additive colored noise. This is a great advantage because the appropriate use 
of the subspace decomposition (possibly in more than the two subspaces as used in this 
thesis), may prove to be an attractive alternative to the existing receiver synchronization 
method in Seaweb. 

B. FUTURE WORK 

There are four major areas that can be identified for future work. First, the 
modeled channel would be more accurate if the statistics used were gathered from the 
real ocean. In this thesis, the variability and spreading imposed upon the impulse 
response were arbitrary. Therefore, the underwater channel may be overestimated or 
underestimated compared to the real ocean. 

Secondly, the algorithms presented here may be implemented in a real 
communication system. The transmitter and receiver groups presented here were 
simplified forms so that the analysis could emphasize the signal processing aspect of the 
scheme. The transmitted signals were considered impulses and the analysis was confined 
to the discrete time domain. A real communication system must account for the effects in 
the frequency domain, as well as in the time domain. 

Thirdly, the loss-of-receiver-synchronization algorithms should be adjusted for 
more eigenvalues in order to accommodate larger shifts. The algorithm presented here 
was proved to equalize adequately when two eigenvalues were used in the simulated 
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doubly spread channel with additive colored noise and a shift of half a bit. More 
eigenvalues should be taken into account if the technique is to replace the existing 
Seaweb synchronization method in a volatile underwater medium. 

Finally, the performance of the proposed receiver structure should be tested at sea. 
Although the system presented here proved to work adequately in a simulated undersea 
environment for various multipath assumptions, noise assumptions, and synchronization 
assumptions, the sea trial is always the ultimate judge of success. 
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APPENDIX A. SOFTWARE USERS MANUAL 


This Appendix serves as a software manual for the doubly spread channel and 
DSSS transceiver structures simulated in this thesis. MATLAB 6.1 release 12 was used to 
implement all aspects of the channel, the transmitter, and the receiver. The MATLAB 
code does not emphasize computational efficiency or real system feasibility but instead 
focuses on the DSP implementation. The choice to use a separate function to initialize the 
modem’s parameters is done to isolate the various parameters from the actual code. It is 
assumed here that the reader has a working knowledge of MATLAB. 

The basic flow of the program will be as illustrated in the Figure 57 below: 



Figure 57. Program Flow Chart 

The impulse response derived from Bellhop is filtered from near zero values and 
saved as ht_stat.mat. In order to inject time variability into the static impulse response the 
m-file time_variable_channel.m must run first. The time-variable channel impulse 
response is saved as ht_var.mat. Subsequently, the time-variable impulse response is 
injected with time spreading using the m-file doubly_spread_channel.m. The doubly 
spread impulse response is saved as ht_gaus.mat. Depending on the case simulated, the 
appropriate impulse response is loaded to the transceiver. 

Once the desired channel is available, the transceiver is setup using the function 


73 








setupjx.m. The function has no inputs or outputs hut saves the transmitting and 
simulation parameters as two individual mat files, which are later loaded by the 
transceiver. These mat files are 

• simulation_params.mat - defines the parameters to be used in the overall 
simulation. These include the number of packets and the values of / Nq 
to be used. 

• transmit_params.mat - defines the parameters needed for the transmitter 
and receiver including bit rate, chip rate, sampling frequency, packet 
length and upsampling ratio. 

The available transceiver functions are the following: 

• tranc _w_noise_only.m - For the simulation and performance testing of the 
receiver structure shown in Figure 26, for the white noise only case. No 
channel input is necessary. 

• tranc_stat_ch_only.m - For the simulation and performance testing of the 
receiver structure illustrated in Figure 30, for multipath only. 

• tranc_stat_ch_ w_noise.m - For the simulation and performance testing of 
the receiver structure represented in Figure 30, for the combined effect of 
multipath and additive white noise. 

• tranc_col_noise_only.m - For the simulation and performance testing of 
the receiver structure represented in Figure 30, adapted for colored noise 
only. No channel input is necessary. 

• tranc_stat_ch_ col_noise.m - For the simulation and performance testing 
of the receiver structure represented in Figure 30, for the combined effect 
of multipath and additive colored noise. 

• tranc_gaus_ch_only.m - For the simulation and performance testing of the 
receiver structure illustrated in Figure 30, for multipath in the doubly 
spread channel. 

• tranc_gaus_ch_ col_noise.m - For the simulation and performance testing 
of the receiver structure represented in Figure 30, for the combined effect 
of multipath in a doubly spread channel and additive colored noise. 
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• tranc_no_sync.m - For the simulation and performance testing of the 
receiver structure represented in Figure 30, adapted for the case that no 
synchronization has been achieved between the receiver and the 
transmitter. This case is only simulated for multipath in a doubly spread 
channel and additive colored noise. 

Each of these transceiver structures saves the simulation results as a separate mat- 
fde. This results mat-fde contains a variable called summary. This variable summarizes 
the simulation errors and simulation parameters. The exact name of the mat-fde 
represents the channel and simulation parameters in the following order: channel name, 
number of packets, bit rate, packet length and channel length. For example, 
color_gaus_channel_l_40_72_ch50.mat means doubly spread channel with additive 
colored noise, 1 packet sent, bit-rate 40 bps, packet length 72 bits and channel length 50 
samples. 

Finally, the function plot_results.m plots the results of the simulation summary as 
BERvs. EjN^ . 
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APPENDIX B. MATLAB CODE 


% ------ -. -- 

% This program introduces time variability on the Bellhop Static impulse response 
% and saves the time variable h(t) as ht var .Run this first or use saved ht var 
% to produce the time variable impulse response. 

%Georgios Pelekanos 
%last revision May 8 2003 

% -------- 

clear 

clc 

%load the time invariable ht from Bellhop 

h=load ('cir_40.mat'); 

tau=h.tau'; 

tau=tau(487:end); 

r=h.r; 

hmat=h.hmat; 

zs=h.zs; 

zr=h.zr; 

tref=h.tref; 

ht_stat=h.hmat(:,501)'; 

ht_stat=ht_stat(487:end);%neglect negative times 
save ht stat ht stat; 

%cleanup values less than 0.03 
for i=l:length(ht_stat); 
if abs(ht_stat(i))<0.03; 

ht_stat(i)=0; 

end 

end 

%find the non-zero locations 
taps = find(ht_stat~=0); 

%neglect the direct path which is stable 
tapsl=taps(9:end); 
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disp('started') 
for kk=l: 1000 
kk 

%Create random lags the length of taps 
lags_aux=round(5 * randn( 1 ,length(taps 1))); 

%The first lag set positive 
if lags_aux(l)<0 

lags_aux( 1 )=abs(lags_aux( 1)); 
end 

%Create time variability as rand lags 
rand_lags=taps 1 +lags_aux; 

%Insert the variability by concatenating the direct and variable multipath 
%portions of taps 
var_taps=[taps(l:8) rand_lags]; 
var_taps=sort(var_taps); 

%Check for Duplicates and re-sort 
for d=l:length(var_taps)-l 
if var_taps(d)==var_taps(d+l) 
var_taps(d)=var_taps(d)-1; 
end 
end 

var_taps=sort(var_taps); 

% %Create uniformly distributed phase 
% phi=unifmd(-pi,pi, 1 ,length(taps)); 

%The new time variable impulse response 

ht_var=zeros(size(ht_stat)); 
for i=l:length(taps) 

ht_var(var_taps(i))=ht_stat(taps(i)); 

end 

%Normalize and truncate to keep the original length 
%ht_var=ht_var/max(max(ht_var)); 
ht_var(kk,:)=[ht_var( 1 :length(ht_stat))]; 
end 
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ht_var=sum(ht_var); 
save ht var ht var; 
disp('end') 

figure; 
subplot(2,l,l) 
plot (tau,ht_var) 

title('Time Variable Impulse Response ') 

xlabel('Time(sec)'); 

ylabel('Magnitude'); 

axis([0 0.085 -0.9 0.3]); 

subplot (2,1,2) 

plot(tau,angle(ht_var)) 

title('Time Variable Impulse Response') 

xlabel('Time(sec)'); 

ylabel('Phase(rad)'); 

axis([0 0.085 -4 4]); 
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% - 

% This program introduces time spreading on the Time Variable impulse response 
% and saves the Gaussian variable h(t) as ht Gaus .Run this second and use the 
% saved ht var to reproduce the time variable impulse response. 

% Georgios Pelekanos 
% Last revision May 8 2003 

% - 


clear 

clc 


% _ 

%load the time invariable impulse response 

% 


h=load('cir_40. mat'); 
tau=h.tau'; 
tau=tau(487:end); 
r=h.r; 

hmat=h.hmat; 

zs=h.zs; 

zr=h.zr; 

tref=h.tref; 

fs=40960; 

ht_stat=h.hmat(:,501)'; 

ht_stat=ht_stat(487:end); 

%cleanup values less than 0.05 
for i=l :length(ht_stat); 
if abs(ht_stat(i))<0.03; 

ht_stat(i)=0; 

end 

end 


% _ 

%load the gold code for spread spectrum 

% _ 

load code; 
tx=code(l:50); 

% tx=sign(randn(l,50)); 

% _ 

%load the time variable impulse response 

% _ 

load ht var; 

%Create a gaussian shaped pulse 
var=[8,16]; 

y=[]; 


figure; 
for i=l:2 

y=exp(-(-30:30).^2/(2*var(i))); 
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%Convolve to shift it to the variable taps position 
ht_gaus=conv(ht_var,y); 

%The first 100-300 samples correspond to the direct path 
%and are set equal to a delta function at the original 
%time variable location 

ht_gaus( 1,100:250)=ht_stat( 1,100:250); 
htgaus (1,600:650)=ht_var( 1,600:650); 

ht_gaus=ht_gaus( 1 :length(ht_var)); 

subplot(2,l,i) 

plot(tau,ht_gaus) 

title(sprintf('Time Domain Representation of Doubly Spread Channel Impulse Response for var= 
%.0f,var(i))) 

xlabel('Time(sec)'); 
ylabel('Magnitude'); 
axis([0 0.015 -0.8 0.3]); 

end 

t = tau; 

figure; 

specgram(ht_stat,256,fs,kaiser(256,5),220) 

title('Spectrogram of a Static Impulse Response') 

xlabel('Time(sec)'); 

ylabel('Frequency(Hz)'); 

colorbar('vert') 

axis([0 0.06 0 15000]); 

figure; 
for i=l:2 


y=exp(-(-30:30).^2/(2*var(i))); 

%Convolve to shift it to the variable taps position 

ht_gaus=conv(ht_var,y); 

ht_gaus( 1,100:250)=ht_stat( 1,100:250); 

ht_gaus(l,600:700)=ht_var(l,600:700); 

%The first 590 samples correspond to the direct path 
%and are set equal to a delta function at the original 
%time variable location 

ht_gaus=ht_gaus( 1 :length(ht_var)); 

subplot(2,l,i) 

specgram(ht_gaus,256,fs,kaiser(256,5),220) 

title(sprintf('Spectrogram of Time Variable and Time Spread Impulse Response for var= 
%.0f,var(i))) 

xlabel('Time(sec)'); 
ylabel('Frequency(Hz)'); 
colorbar('vert') 
axis([0 0.06 0 15000]); 
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hold on; 
end 


yx=conv(tx,ht_stat); 
yx_var=conv(tx,ht_var); 
yx_gaus=conv(tx,ht_gaus); 


HT=fft(ht_stat); 

HT_var=fft(ht_var); 

HT_gaus=fft(ht_gaus); 

F= ( 0: (fs/2)/4090: (fs/2)*4089/4090 ); 

figure; 

subplot(2,l,l) 

plot(F,20*logl0(fftshift(abs(HT(l:4090))))) 
title('Spectral Content of a Static Impulse Response ') 
xlabel('Frequency(Hz)'); 
ylabel('Magnitude(dB)'); 

subplot(2,l,2) 

plot(F,20*logl0(fftshift(abs(HT_gaus(l:4090))))) 

title(sprintf('Spectral Content of a Doubly Spread Impulse Response for var= %.0f,var(2))) 

xlabel('Frequency(Hz)'); 

ylabel('Magnitude(dB)'); 


Yx=fft(yx); 

Y x=fft(yx_var); 

Y x_gaus=fft(yx_gaus); 

F= ( 0: (fs/2)/4090: (fs/2)*4089/4090 ); 

figure; 

subplot(2,l,l) 

plot(F,20* log 10(fftshift(abs(Y x(1:4090))))) 
title('Spread Spectmm over a Static Channel') 
xlabel('Frequency(Hz)'); 
ylabel('Magnitude(dB)'); 

subplot(2,l,2) 

plot(F,20*logl0(fftshift(abs(Yx_gaus(l:4090))))) 

title(sprintf('Spread Spectmm over a Doubly Spread Impulse Response for var= %.0f ,var(2))) 

xlabel('Frequency(Hz)'); 

ylabel('Magnitude(dB)'); 

stat_corr=xcorr(ht_stat); 

gaus_corr=xcorr(ht_gaus); 


figure; 

subplot(2,l,l) 

plot (stat corr) 

axis ([1000 7000-0.1 1]) 

subplot(2,l,2) 
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plot (gaus corr) 

axis ([1000 7000 -0.3 2.5]) 

save ht gaus ht gaus; 

% [Pxx_stat,wl] = pburg(yx,12,512,40960,'onesided'); 

% [Pxx_gaus,w2] = pburg(yx_gaus,12,512,40960,'onesided'); 

[Pxx_stat,wl] = pwelch(ht_stat,[],[],40960) 

[Pxx_gaus,w2] = pwelch(ht_gaus,[],[],40960) 

figure; 

subplot(2,l,l) 

psdplot(Pxx_stat,wl ,'Hz') 

subplot(2,l,2) 

psdplot(Pxx_gaus,w2,'Hz') 
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function setup tx 

clear 

clc 

% 

% This function sets up all transmit parameters 
% and saves them as separate *.mat files 
% Two separate mat files are generated 
% 1. simulation_params.mat 
% 2. transmit_params. mat 
% 

% developed by Georgios Pelekanos February 2003 
% last modified 27/2 

%Transmit parameters 
Rb = 40; 

Rc = 5000; % chip rate 

Rs = 40960; % sampling rate 

fs = Rs; % sampling frequency 

Tb = 1/Rb; % bit length 

Tc = 1/Rc; % chip length 

Ts = 1/fs; % sampling interval 

fcarrier = 12000; % carrier frequency 

samples_per_bit = Rs/Rb; % samples per bit 

samples_per_chip = Rs/Rc; % samples per chip 

chips_per_bit = Rc/Rb; % chips per bit 

packet length = 72; % number of dbits transmitted per packet 

P=Tb/Tc; % Upsampling required to match chip sequence 

% save settings to "transmit_params.mat" file 

save transmit_params Rb Rc Rs fs Tb Tc Ts fcarrier P ... 

samples_per_bit samples_per_chip chips_per_bit... 
packet length ; 


% Channel Simulation Parameters 

N = 1; %number of packets 

EbNo_dB = 35; %EbNo values to use 

% save channel settings to "simulation_params.mat" file 

save simulation_params N EbNo dB; 
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% Tranceiver v. 1 
% WHITE NOISE ONLY 

% developed by Georgios Pelekanos February 2003 
% last modified 15/4 

clear; 

clc; 

tic; % starts clock to measure simulation run time 

disp('started'); 

load code; % load PN chipping sequence code.mat 

load transmit_params; % load transmit parameters 
load simulation_params; % load simulation parameters 
load htstat; 

EbNo_dB=linspace(l , EbNo_dB, 10); 

for pp=l: 1000 
PP 

for kl = l:length(EbNo_dB); %loop through each of the EbNo values 


% 


%%%%%% %%%%% TRANSMITTER SECTION %%%%%%%%%%%%%%%%% 
% 


% - 

% Generate a differential bit sequence of 1,-1 of 
% length=N*packet_length simulating the signal 
%------ 

signal_an=sign(randn( 1 ,N*packet_length)); 

% .. —- . —- -- 

% Upsample by P=Tb/Tc 
%--- 


upsamples =[l,ones(l,P-l)]; 

for i=l:N*packet_length 

matrix):,i) = signal_an(i).*upsamples' ; 
end 

[m,n]=size(matrix); 

up_sampled_an=reshape(matrix, 1 ,m*n); 


% -- . .. - . - -- 

% Chip the data using as c[n] a multiple/portion of original gold code 
% first P samples of gold code only) 

% -- . - ----- -- 


% Determine repetitions of orignal gold code required 
% to chip entire upsampled series 

totalchipsneeded = length(signal_an).*chips_per_bit; 
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numberofrepetitions = ceil(total_chips_needed./P); 

c_n = repmat(code(l:P), 1, number_of_repetitions);%long PN sequence 

%produce the transmit signal by multiplication 
%which is the same as c(n)*up_sampled_an 

tx_signal=c_n. *up_sampled_an; 


% ------- 

% Calculate the signal power and associated Eb/No 
%-----.. 


tx_signal_power = (sum(tx_signal.^2))./length(tx_signal); 

% generate the noise power associated with EbNo 

EbNo(kl) = 10."'(EbNo_dB(kl)/10); 

Eb = tx_signal_power.*Tb; 

No = Eb./EbNo(kl); 
sigma_n = sqrt(No./(2*Ts)); 
noise_power = sigma_n.^2; 

SNRdB = 10.*logl0(tx_signal_power./noise_power); 


% 



CHANNEL SECTION %%%%%%%%%%%% 


% ---- 

% Add Gaussian White noise no multipath channel 
%.-.—-.-.. 


noisestd = sqrt(noise_power); 

y_n = txsignal + randn(l,length(tx_signal))*noisestd; 


% 


%%%%%%%%%%%%%% RECEIVER SECTION 
% 


%%%%%%%%%%%%%%%% 


% Create a replica of the transmitted code such that c(n)*c(-n)=delta(n) 

% determine repetitions/portion of orignal gold code required 
% to decode entire output 

numberofrepetitions = ceil(length(y_n)/P); 

c minus n = repmat(code(l:P), 1, number_of_repetitions);%long PN sequence 

%Decode the sequence using the replica and assuming synchronization 

decode_rx= y_n( 1: m*n). *c_minus_n( 1 :m* n) ;%decode 

%The channel estimation algorithm 

% Delay by z^-1 Downsample by P new revised 3/4 works 
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decode_rx=reshape(decode_rx,P,length(decode_rx)/P); 

%Use maximization matched filter for white noise 

R_y=decode_rx* decoderx'; 

[Eig_Vectors,Eig_Values] = eig(R_y); 
lambda = diag(Eig_Values); 
f_bar=Eig_Vectors(:,P); 

% lambda=-l/norm(code(l:P)); 

% f_bar=lambda*code(l:P)'; 


% --- 

%Apply the matched filter 
%--- 


a_hat_n=f_bar'*decode_rx; 
a_hat_n=sign(a_hat_n); 


% --- 

%Calculate Errors 

% .-. 

counter=find(signal_an-a_hat_n~=0); 

if length(counter)<= 10^-6 
channel_errors(kl)= 10^-6; 
else 

channel_errors(kl)=length(counter); 

end 

summary(pp,:)=channel_errors./packet_length; 

end 

errors=mean( summary); 

overall=[EbNo_dB;errors]; 

% save the simulation results 

save montecarlo_white_only_144_40_50 overall; 

disp ('finished'); 
t = toe 


87 







% Tranceiver v.2 FINAL 
% MULTIPATH STATIC CHANNEL ONLY 
% developed by Georgios Pelekanos February 2003 
% last modified 15/4 

clear; 

clc; 

tic; % starts clock to measure simulation run time 

disp('started'); 

load code; % load PN chipping sequence code.mat 

load transmit_params; % load transmit parameters 
load simulation_params; % load simulation parameters 
load htstat; 


for kl = 1:100; %Repeat the simulation 50 times 


% 


%%%%%%%%%%%%%% TRANSMITTER SECTION %%%%%%%%%%%%%%%%% 
% 


% . —- . - -- 

% Generate a differential bit sequence of 1,-1 of 
% length=N*packet_length simulating the signal 
%.—-.... 


signal_an=sign(randn( 1 ,N*packet_length)); 


%—-.—-.-. 

% Upsample by P=Tb/Tc 
%.... 


upsamples =[l,ones(l,P-l)]; 

for i=l:N*packet_length 

matrix(:,i) = signal_an(i).*upsamples' ; 
end 

[m,n]=size(matrix); 

up_sampled_an=reshape(matrix, 1 ,m*n); 


% --- 

% Chip the data using as c[n] a multiple/portion of original gold code 
% first P samples of gold code only) 

% --- 


% Determine repetitions of orignal gold code required 
% to chip entire upsampled series 

totalchipsneeded = length(signal_an).*chips_per_bit; 
numberofrepetitions = ceil(total_chips_needed./P); 
c_n = repmat(code(l:P), 1, number of repetitions); 
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%produce the tranmit signal by multiplication 
%which is the same as c(n)*up_sampled_an 

tx_signal=c_n. *up_sampled_an; 


% 


%%%%%%%%%%%%%%%% 


% 


CHANNEL SECTION %%%%%%%%%%%%%%% 


%- ... .- . - . . 

%convolve with the channel to produce y(n) 
%----- 


y_n = conv(tx_signal,ht_stat(590:end)); 


% 


%%%%%%%%%%%%%%%%% RECEIVER SECTION %%%%%%%%%%%%%%%%%% 


% 


% Create a replica of the transmitted code such that c(n)*c(-n)=delta(n) 

% determine repetitions/portion of orignal gold code required 
% to decode entire output 

numberofrepetitions = ceil(length(y_n)/P); 

c minus n = repmat(code(l:P), 1, number_of_repetitions);%long PN sequence 


%Decode the sequence using the replica and assuming S5mchronization 

decode_rx= y_n( 1: m*n). *c_minus_n( 1 :m* n) ;%decode 

%The channel estimation algorithm 

% Delay by z^-1 Downsample by P new revised 3/4 works 

decode_rx=reshape(decode_rx,P,length(decode_rx)/P); 

%Use maximization matched filter for other than w.noise 

R_y=decode_rx* decoderx'; 

[Eig_Vectors,Eig_Values] = eig(R_y); 
lambda = diag(Eig_Values); 
f_bar=Eig_Vectors(:,P); 


% --- 

%Apply the matched filter 
%...--- 


a_hat_n=f_bar'*decode_rx; 
a_hat_n=sign(a_hat_n); 
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% --- 

%Calculate Errors 

% -- 

counter=find(signal_an-a_hat_n~=0); 
if length(counter)<= 10^-6 
channel_errors(kl)= 10^-6; 
else 

channel_errors(kl)=length(counter); 

end 

percent_ch_errors(kl)=channel_errors(kl)/ (N*packet_length)* 100; 

% generate a summary of all the simulation results 

summary(:,kl) = [N*packet_length, kl, channel_errors(kl),percent_ch_errors(kl)]'; 
end 

% save the simulation results 
save simulation static only summary 

disp ('finished'); 
t = toe 
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% Tranceiver v.3 FINAL 
% MULTIPATH + WHITE NOISE 
% developed by Georgios Pelekanos February 2003 
% last modified 15/4 

clear; 

clc; 

tic; % starts clock to measure simulation run time 

disp('started'); 

load code; % load PN chipping sequence code.mat 

load transmit_params; % load transmit parameters 
load simulation_params; % load simulation parameters 
load htstat; 

EbNo_dB=linspace(l , EbNo_dB, 10); 

for pp= 1:1000 
PP 

for kl = l:length(EbNo_dB); %loop through each of the EbNo values 


% 


%%%%%%%%%%%%%%%% TRANSMITTER SECTION %%%%%%%%%%%%%%%% 
% 


% - 

% Generate a differential bit sequence of 1,-1 of 
% length=N*packet_length simulating the signal 
%---- 

signal_an=sign(randn( 1 ,N*packet_length)); 

% .. —- . —- -- 

% Upsample by P=Tb/Tc 
%--- 


upsamples =[l,ones(l,P-l)]; 

for i=l:N*packet_length 

matrix):,i) = signal_an(i).*upsamples' ; 
end 

[m,n]=size(matrix); 

up_sampled_an=reshape(matrix, 1 ,m*n); 


% -- . .. - . - -- 

% Chip the data using as c[n] a multiple/portion of original gold code 
% first P samples of gold code only) 

% ..- . - ----- - 


% Determine repetitions of orignal gold code required 
% to chip entire upsampled series 

totalchipsneeded = length(signal_an).*chips_per_bit; 
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numberofrepetitions = ceil(total_chips_needed./P); 

c_n = repmat(code(l:P), 1, number_of_repetitions);%long PN sequence 

%produce the tranmit signal by multiplication 
%which is the same as c(n)*up_sampled_an 

tx_signal=c_n. *up_sampled_an; 


% ------- 

% Calculate the signal power and associated Eb/No 
%----- 


tx_signal_power = (sum(tx_signal.^2))./length(tx_signal); 

% generate the noise power associated with EbNo 

EbNo(kl)= 10."'(EbNo_dB(kl)/10); 

Eb = tx_signal_power.*Tb; 

No = Eb./EbNo(kl); 
sigma_n = sqrt(No./(2*Ts)); 
noise_power = sigma_n.^2; 

SNRdB = 10.*logl0(tx_signal_power./noise_power); 


% 


%%%%%%%%%%%%%%%% CHANNEL SECTION %%%%%%%%%%%%%%%%%%% 


%—------- 

%convolve with the channel to produce y(n) 
%------ 

y_n = conv(tx_signal,ht_stat(590:640)); 

% . - . - .. 

% Add Gaussian White Noise 
%------- 

noisestd = sqrt(noise_power); 

y_n = y_n + randn(l,length(y_n))*noisestd; 

% 


%%%%%%%%%%%%%%%% RECEIVER SECTION %%%%%%%%%%%%%%%%%%% 


% 


% Create a replica of the transmitted code such that c(n)*c(-n)=delta(n) 

% determine repetitions/portion of orignal gold code required 
% to decode entire output 

numberofrepetitions = ceil(length(y_n)/P); 

c minus n = repmat(code(l:P), 1, number_of_repetitions);%long PN sequence 
%Decode the sequence using the replica and assuming synchronization 
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decode_rx= y_n( 1: m*n). *c_minus_n( 1 :m* n) ;%decode 
%The channel estimation algorithm 
% Delay by z^-1 Downsample by P new revised 3/4 works 
decode_rx=reshape(decode_rx,P,length(decode_rx)/P); 
%Matched filter SNR maximization Last revised 17/3 works 
if channeltype == 0; 

%Use maximization matched filter for white noise 

lambda=- l/norm(c_minus_n( 1 :P)); 
f_bar=lambda*c_minus_n(l :P)'; 

else 

%Use maximization matched filter for other than w.noise 

R_y=decode_rx*decode_rx'; 

[Eig_Vectors,Eig_Values] = eig(R_y); 
lambda = diag(Eig_Values); 
f_bar=Eig_Vectors (: ,P); 

end 

%—- . - . - . 

%Apply the matched filter 
%.... 


a_hat_n=f_bar'*decode_rx; 
a_hat_n=sign(a_hat_n); 


% -- 

%Calculate Errors 

% --- 

counter=find(signal_an-a_hat_n~=0); 
if length(counter)<= 10^-6 
channel_errors(kl)= 10^-6; 
else 

channel_errors(kl)=length(counter); 

end 

summary(pp,:)=channel_errors./packet_length; 

end 

errors=mean( summary); 
overall=[EbNo_dB;errors]; 

% save the simulation results 
save montecarlo_white_only_144_40_50 overall; 
disp ('finished'); 
t = toe 
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% Tranceiver v.4 
% COLORED NOISE ONLY 
% developed by Georgios Pelekanos February 2003 
% last modified 5/12 

clear; 

clc; 

tic; % starts clock to measure simulation run time 

disp('started'); 

load code; % load PN chipping sequence code.mat 

load transmit_params; % load transmit parameters 
load simulation_params; % load simulation parameters 


% 


%%%%%%%%%%%%% TRANSMITTER SECTION %%%%%%%%%%%%%%%%%%% 
% 


EbNo_dB=linspace(l , EbNo_dB, 10); 

for pp= 1:1000 
PP 

for kl = l:length(EbNo_dB); %loop through each of the EbNo values 

% --- 

% Generate a differential bit sequence of 1,-1 of 
% length=N*packet_length simulating the signal 
%--- 

signal_an=sign(randn( 1 ,N*packet_length)); 

% . - . - . - . —- 

% Upsample by P=Tb/Tc 
%----- 

upsamples =[l,ones(l,P-l)]; 

for i=l:N*packet_length 

matrix):,i) = signal_an(i).*upsamples' ; 
end 

[m,n]=size(matrix); 

up_sampled_an=reshape(matrix, 1 ,m*n); 

%—-.-.— 

% Chip the data using as c[n] a multiple/portion of original gold code 
% first P samples of gold code only) 

% ......... — 


% Determine repetitions of orignal gold code required 
% to chip entire upsampled series 

totalchipsneeded = length(signal_an).*chips_per_bit; 

numberofrepetitions = ceil(total_chips_needed./P); 

c_n = repmat(code(l:P), 1, number_of_repetitions);%long PN sequence 
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%produce the tranmit signal by multiplication 
%which is the same as c(n)*up_sampled_an 

tx_signal=c_n. *up_sampled_an; 

% ----- 

% Calculate the signal power and associated Eb/No 
%—-----.-.- 


tx_signal_power = (sum(tx_signal.^2))./length(tx_signal); 

% generate the noise power associated with EbNo 

EbNo(kl) = 10.^(EbNo_dB(kl)/10); 

Eb = tx_signal_power.*Tb; 

No = Eb./EbNo(kl); 
sigma_n = sqrt(No./(2*Ts)); 
noise_power = sigma_n.^2; 

SNRdB = 10.*logl0(tx_signal_power./noise_power); 

% 


%%%%%%%%%%%%%%%% CHANNEL SECTION %%%%%%%%%%%%%%%%%%% 


% 


%—-.—-.- 

% Add Colored noise 
% -.—- 


noisestd = sqrt(noise_power); 
white_noise=randn( 1 ,length(tx_signal)); 
color_noise=filter(l,[l -0.9 0.87 -0.89 0.9 -O.88],white_noise); 
y_n = tx signal + color_noise*noisestd ; 


% 


%%%%%%%%%%%%%%% RECEIVER SECTION %%%%%%%%%%%%%%%%%%% 


% 


% Create a replica of the transmitted code such that c(n)*c(-n)=delta(n) 

% determine repetitions/portion of orignal gold code required 
% to decode entire output 

numberofrepetitions = ceil(length(y_n)/P); 

c minus n = repmat(code(l:P), 1, number_of_repetitions);%long PN sequence 

%Decode the sequence using the replica and assuming synchronization 

decode_rx= y_n(l :m*n).*c_minus_n(l :m*n);%decode 

%The channel estimation algorithm 

% Delay by z^-1 Downsample by P new revised 3/4 works 

decode_rx=reshape(decode_rx,P,length(decode_rx)/P); 
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%Use maximization matched filter for colored noise 
R_y=decode_rx* decoderx'; 

color_noise=[0,color_noise]; 

R_n=color_noise*color_noise'; 

[Eign, Vain] =eig(R_n); 
Rnsq_inv=inv(sqrtm(Val_n)); 

R_n_sq_inv=Eig_n * Rnsq_inv* Eign'; 
R_mah=R_n_sq_inv* R_y * R_n_sq_inv; 

[Eigmah, Valmah] =eig(R_mah); 
[Eig_mah,Val_mah] = sortem(Eig_mah,Val_mah); 
h_pr=Eig_mah(:, 1); 
h_opt=R_n_sq_inv*h_pr; 

if h_opt(2:4)<0 
h_opt=-h_opt; 
end 


% - 

%Apply the matched filter 
%- 


a_hat_n=h_opt'*decode_rx; 
a_hat_n=sign(a_hat_n); 


% .-. 

%Calculate Errors 

% ---— 

counter=find(signal_an-a_hat_n~=0); 
if length(counter)<= 10^-6 
channel_errors(kl)= 10^-6; 
else 

channel_errors(kl)=length(counter); 

end 

summary(pp,:)=channel_errors./packet_length; 

end 

errors=mean( summary); 

overall=[EbNo_dB;errors]; 

% save the simulation results 

save montecarlo_color_only_144_40_50 overall; 
disp ('finished'); 
t = toe 
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% Tranceiver v.5 

% COLORED NOISE and STATIC CHANNEL or DOUBLY SPREAD CHANNEL 
% for static load ht stat 
% for doubly spread load ht gaus 
% developed by Georgios Pelekanos Eebruary 2003 
% last modified 5/12 

clear; 

clc; 

tic; % starts clock to measure simulation run time 

disp('started'); 

load code; % load PN chipping sequence code.mat 

load transmit_params; % load transmit parameters 
load simulation_params; % load simulation parameters 
load ht gaus; 

% 


%%%%%%%%%%%%% TRANSMITTER SECTION %%%%%%%%%%%%%%%%%%% 
% 


EbNo_dB=linspace(l , EbNo_dB, 8); 

for pp=l: 1000 
PP 

for kl = 1 :length(EbNo_dB); %loop through each of the EbNo values 


% . - --- - -- 

% Generate a differential bit sequence of 1,-1 of 
% length=N*packet_length simulating the signal 
%----- 


signal_an=sign(randn( 1 ,N*packet_length)); 


% --- 

% Upsample by P=Tb/Tc 
%- 


upsamples =[l,ones(l,P-l)]; 

for i=l:N*packet_length 

matrix(:,i) = signal_an(i).*upsamples' ; 
end 

[m,n]=size(matrix); 

up_sampled_an=reshape(matrix, 1 ,m*n); 


% .-- ------ - - 

% Chip the data using as c[n] a multiple/portion of original gold code 
% first P samples of gold code only) 

% --- 
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% Determine repetitions of orignal gold code required 
% to chip entire upsampled series 

totalchipsneeded = length(signal_an).*chips_per_bit; 

numberofrepetitions = ceil(total_chips_needed./P); 

c_n = repmat(code(l:P), 1, number_of_repetitions);%long PN sequence 


%produce the tranmit signal by multiplication 
%which is the same as c(n)*up_sampled_an 

tx_signal=c_n. *up_sampled_an; 

% - - — --- - 

% Calculate the signal power and associated Eb/No 
%—-.—-.-.-.-.- 


tx_signal_power = (sum(tx_signal.^2))./length(tx_signal); 

% generate the noise power associated with EbNo 

EbNo(kl) = 10."'(EbNo_dB(kl)/10); 

Eb = tx_signal_power.*Tb; 

No = Eb./EbNo(kl); 
sigma_n = sqrt(No./(2*Ts)); 
noise_power = sigma_n.^2; 

SNRdB = 10.*logl0(tx_signal_power./noise_power); 

% 


%%%%%%%%%%%%%%CHANNEL SECTION %%%%%%%%%%%%%%%%%%%%% 


% 


% --- 

%convolve with the channel to produce y(n) 
%—------- 

% y_n = conv(tx_signal,ht_stat(590:640)); 
y_n = conv(tx_signal,ht_gaus(90:190)); 

% .. .- .. . . 

% Add Colored Noise 

% ------ 


noisestd = sqrt(noise_power); 
white_noise=randn( 1 ,length(y_n)); 
color_noise=filter(l,[l -0.9 0.4 -0.3 0.2 0.1],white_noise); 

y_n = y_n + color_noise*0.7*noisestd ; 


% 


%%%%%%%%%%%%%%% RECEIVER SECTION %%%%%%%%%%%%%%%%%%%% 


% 
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% Create a replica of the transmitted code such that c(n)*c(-n)=delta(n) 
% determine repetitions/portion of orignal gold code required 
% to decode entire output 


numberofrepetitions = ceil(length(y_n)/P); 

c minus n = repmat(code(l:P), 1, number_of_repetitions);%long PN sequence 

%Decode the sequence using the replica and assuming synchronization 

decode_rx= y_n( 1: m*n). *c_minus_n( 1 :m* n) ;%decode 
%The channel estimation algorithm 
% Delay by z^-1 Downsample by P new revised 3/4 works 
decode_rx=reshape(decode_rx,P,length(decode_rx)/P); 


%Use maximization matched filter for colored noise 
R_y=decode_rx* decoderx'; 

color_noise=[0,color_noise]; 

R_n=color_noise*color_noise'; 

[Eign, Vain] =eig(R_n); 
Rnsq_inv=inv(sqrtm(Val_n)); 

R_n_sq_inv=Eig_n * Rnsq_in v* Eign'; 
R_mah=R_n_sq_inv* R_y * R_n_sq_inv; 

[Eigmah, Valmah] =eig(R_mah); 
[Eig_mah,Val_mah] = sortem(Eig_mah,Val_mah); 
h_pr=Eig_mah(:, 1); 
h_opt=R_n_sq_inv*h_pr; 

if h_opt(2:4)<0 
h_opt=-h_opt; 
end 


% ---- 

%Apply the matched filter 
%--- 


a_hat_n=h_opt'*decode_rx; 
a_hat_n=sign(a_hat_n); 


% - 

%Calculate Errors 

% . 

counter=find(signal_an-a_hat_n~=0); 

if length(counter)<= 10^-6 
channel_errors(kl)= 10^-6; 
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else 

channel_errors(kl)=length(counter); 

end 

summary(pp,:)=channel_errors./packet_length; 

end 

errors=niean( summary); 
overall=[EbNo_dB;errors]; 

% save the simulation results 

save montecarlo_gaus_color_1242_40_100 overall; 


% %-.-.—-- 

% % Plot overall simulation results 



% es = summary(5,:)./(summary(2,:)); 


for i=l:length(EbNo_dB) 
if overall(2,i)==0 
overall(2,i)=10^-6 
end 
end 

figure; 

semilogy(overall(l,:),overall(2,:),'-','MarkerSize',8); 

grid 

ylim([10"'-5 l]);xlim([l 22]); 
ylabel('Bit Error Rate','fontsize',13); 
xlabel('Eb/No (dB)','fontsize', 13); 
set(gca,'fontsize', 13); 

disp ('finished'); 
t = toe 
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% Blind Tranceiver when no sync 
%Tranceiver No Synch with MonteCarlo 
% developed by Georgios Pelekanos April 2003 
% last modified 5/12 

clear; 

clc; 

tic; % starts clock to measure simulation run time 

disp('started'); 

load code; % load PN chipping sequence code.mat 

load transmit_params; % load transmit parameters 
load simulation_params; % load simulation parameters 
load htstat; 
load htgaus; 

EbNo_dB=linspace(l ,EbNo_dB,8); 

for pp=l: 1000 
PP 

for kl = l:length(EbNo_dB); %loop through each of the EbNo values 
% 


%%%%%%%%%%%%%% TRANSMITTER SECTION %%%%%%%%%%%%%%%%%% 
% 


% ----- —- 

% Generate a differential bit sequence of 1,-1 of 
% length=N*packet_length simulating the signal 
%.-.-. 

signal_an=sign(randn( 1 ,N*packet_length)); 

% . - .. .—- 

% Upsample by P=Tb/Tc 
%--- 


upsamples =[l,ones(l,P-l)]; 

for i=l:N*packet_length 

matrix(:,i) = signal_an(i).*upsamples' ; 
end 

[m,n]=size(matrix); 

up_sampled_an=reshape(matrix, 1 ,m*n); 


%—-.- . — 

% Chip the data using as c[n] a multiple/portion of original gold code 
% first P samples of gold code only) 

% . - .. . . . . 


% Determine repetitions of orignal gold code required 
% to chip entire upsampled series 

totalchipsneeded = length(signal_an).*chips_per_bit; 

numberofrepetitions = ceil(total_chips_needed./P); 

c_n = repmat(code(l:P), 1, number_of_repetitions);%long PN sequence 
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%produce the tranmit signal by multiplication 
%which is the same as c(n)*up_sampled_an 

tx_signal=c_n. *up_sampled_an; 

% . - . —- . 

% Calculate the signal power and associated Eb/No 

% --- 

tx_signal_power = (sum(tx_signal.^2))./length(tx_signal); 

% generate the noise power associated with EbNo 

EbNo(kl) = 10."'(EbNo_dB(kl)/10); 

Eb = tx_signal_power.*Tb; 

No = Eb./EbNo(kl); 
sigma_n = sqrt(No./(2*Ts)); 
noise_power = sigma_n.^2; 

SNRdB = 10.*logl0(tx_signal_power./noise_power); 


% 


%%%%%%%%%%%%%%%% CHANNEL SECTION %%%%%%%%%%%%%%%%%%% 


% 


% g=ht_stat(590:640); 
g=ht_gaus(90:140); 

% -- - - - --- - 

%convolve with the channel to produce y(n) 

% - 

y_n = conv(tx_signal,g); 

% - 

% Add Gaussian White Noise 
%--- 

noisestd = sqrt(noise_power); 

y_n = y_n + randn(l,length(y_n))*0.7*noisestd; 

% 


%%%%%%%%%%%%%%%% RECEIVER SECTION %%%%%%%%%%%%%%%%%%% 


% 


% To simulate lack of synchronization we add random initial time tO 
% uniformly distributed between 0 and P. 

tO=round(rand*P/2)+1; 
y_n=y_n(t0: length(y_n)); 

% Create a replica of the transmitted code such that c(n)*c(-n)=delta(n) 
% determine repetitions/portion of orignal gold code required 
% to decode entire output 
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numberofrepetitions = ceil(length(y_n)/P); 

c minus n = repmat(code(l:P), 1, number_of_repetitions);%long PN sequence 

%Decode the sequence using the replica and assuming synchronization 

decode_rx= y_n( 1: m*n). *c_minus_n( 1 :m* n) ;%decode 

%The channel estimation algorithm 

% Delay by z^-1 Downsample by P new revised 3/4 works 

decode_rx=reshape(decode_rx,P,length(decode_rx)/P); 

%Use maximization matched filter for other than w.noise 

R_y=decode_rx* decoderx'; 

[Eig_Vectors,Eig_Values] = eig(R_y); 

Eig_Values=diag(Eig_Values); 

if EigV alues(P-l)>3/10* EigV alues (P) 

% Assume there is no sync 
e=Eig_V ectors(:,P-l:P)'* decoderx; 

E=[e(:, 1:length(e)-1); e(: ,2:length(e))]; 

Re=E*E'; 

[Ve, De]=eig(Re); 
lam_e=diag(De); 

I=find(lam_e==min(abs(lam_e))); 

me=Ve(:,I); 

ml=me(l:2); 

m2=-me(3:4); 


% --- 

% Estimate the sequence 

% . - . - . 

a_hat_n=[ml ',m2'] *E; 
a_hat_n=sign(a_hat_n); 

else 

% Assume there is sync 
f_bar=Eig_Vectors(:,P); 

0 /„- 

%Apply the matched filter 

% --- 

a_hat_n=f_bar'*decode_rx; 

a_hat_n=sign(a_hat_n); 

end 


% ...- 

%Calculate Errors 

% --- 

diff=length(a_hat_n)-length(signal_an); 

ifdiff>0 

counter=find(signal_an-a_hat_n(diff:length(signal_an)+diff-l)~=0); 
elseif diff<0 
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diff=abs(diff); 

counter=find(signal_an(diff:length(a_hat_n)+diff-l)-a_hat_n~=0) 

else 

counter=find(signal_an-a_hat_n~=0); 

end 

if length(counter)<= 10^-6 
channel_errors(kl)= 10^-6; 
else 

channel_errors(kl)=length(counter); 

end 

sumniary(pp,:)=channel_errors./packet_length; 

end 

errors=niean( summary); 

overall=[EbNo_dB;errors]; 

% save the simulation results 

save montecarlo_no_sync_144_40 overall; 

disp ('finished'); 
t = toe 
break 
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fiinction plotresult 


clear 

clc 

% 

% This function sets plots the results 

% 

% developed by Georgios Pelekanos February 2003 
% last modified 27/4 


load montecarlo_no_sync_144_40 ; 


% % -- 

% % Plot overall simulation results 

% % ....- 

% 


for i=l :length(EbNo_dB) 
if overall(2,i)==0 
overall(2,i)=10^-6 
end 
end 

figure; 

semilogy(overall(l,:),overall(2,:),'-','MarkerSize',8); 

grid 

ylim([10^-5 1]); xlim([l 23]); 
ylabel('Bit Error Rate','fontsize',13); 
xlabel('Eb/No (dB)','fontsize', 13); 
set(gca,'fontsize', 13); 
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